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1.0  INTRODUCTION 

The  design  criteria  for  high-speed  ship  propellers  involve 
tradeoffs  between  efficiency  and  cavitation,  and  strength  and 
vibration  of  the  propeller.  Operating  in  ship  wakes  at 
relatively  low  cavitation  numbers,  the  propeller  will,  at  least 
intermittently,  cavitate,  inducing  erosion  of  the  blades,  loss 
of  efficiency,  noise,  vibration,  and  occasionally  structural 
failure  of  nearby  plating.  These  harmful  effects  are  mainly  due 
to  the  collapse  of  unsteady  cavities.  These  include  individual 
bubbles  as  well  as  sheet  cavities  and  "clouds"  (Figure  1.1). 

Adequate  and  increasingly  sophisticated  theories  for 
individual  bubble  growth  and  collapse  exist  (see  the  reviews  by 
Plesset  and  Prosperetti,  1977,  and  Hammitt,  1980).  While  the 
transition  to  sheet  cavity  is  not  well  understood,  a  large 
number  of  experimental  observations  of  sheet  cavitation  are 
available  (Shen  and  Peterson,  1979,  Bark  and  Barlekom,  1979,  and 
Bark,  1985).  A  steady,  then  unsteady,  theory  for  sheet 
cavitation  was  recently  developed  (Tulin,  1980,  and  Tulin  and 
Hsu,  1980).  Downstream  of  a  "steady"  sheet  cavity  a  region  of 
high  population  of  tiny  bubbles  can  be  observed  and  is 
especially  known  to  be  associated  with  erosion.  However,  as 
concluded  from  observations  by  Tanabayachi  and  Chiba  (1977),  an 
unsteady  sheet  cavity  is  required  for  the  formation  of  coherent 
clouds  of  very  fine  bubbles.  These  clouds  are  either  detached 
from  the  frothy  mixture  at  the  trailing  end  of  the  unsteady 
sheet,  or  generated  in  a  finite  region  of  the  liquid  downstream 
of  the  unsteady  sheet  where  significant  fluctuating  pressures 
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As  the  pressures  generated  by  single  bubble  collapse  are 
not  strong  enough  to  explain  the  intense  erosion  in  the  subject 
region  and  the  high  forces  needed,  for  example,  to  bend  the 
trailing  edge,  cloud  cavitation  has  been  held  responsible  since 
Van  Manen's  (1963)  work.  This  is  supported  experimentally  by  a 
very  close  correlation  between  the  dynamics  of  these  clouds  and 
the  sharpest  and  highest  pressure  pulses  detected  on  an 
oscillating  hydrofoil  (Bark  and  Barlekom,  1979).  Similar 
phenomena  have  been  observed  with  ultrasonic  cavitation  (Hanson 
and  M^rch,  1 980  )  . 

Apart  from  some  information  on  the  frequency  of  generation 
of  cloud  cavitation,  the  experimental  observations  and 
measurements  are  very  qualitative  and  do  not  allow  at  the 
present  time  any  quantitative  predictions.  In  addition,  the 
lack  of  understanding  of  the  dynamics  of  such  cavities  makes  it 
impossible  to  explain  any  scaling  effects  and  to  correct  for 
them.  Theoretical  and  fundamental  studies  are  thus  needed  as 
guidance  for  future  design  and  experimentation. 

To  our  knowledge,  since  the  early  work  of  Van  Wijngaarden 
(1964)  only  a  few  publications  by  M^rch  (1977,  1980,  and  1982) 
and  Hanson  and  M^rch  (1980)  have  dealt  theoretically  with  the 
problem  of  "collective  bubbles  collapse"  or  "cavity  cluster 
collapse."  However,  a  large  amount  of  literature  has  been 
devoted  to  the  modeling  of  bubble-liquid  mixture  behavior  using 
either  a  continuum  medium  approach  or  a  "two-fluid"  approach 
(Zwick,  1959;  Van  Wijngaarden,  1972,  1976,  1980,  and  1982; 

Zuber,  1964;  Ishii,  1975).  This  subject  has  recently  regained  a 
lot  of  interest  as  is  shown  by  several  very  recent  publications 
(D'Agostino  and  Brennen;  1983,  Rubinstein,  1985;  Caflish,  et. 
al,  1985a-b;  Ng  and  Ting,  1986;  Miksis  and  Ting,  1 986; . . .etc . ) 
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In  order  to  explain  the  phenomenon  of  propeller  blades  bent  at 
the  trailing  edge.  Van  Wijngaarden  (1964)  considered  the  case  of 
a  uniform  layer  of  cavities  on  a  solid  wall.  He  studied  its 
unidimensional  collective  collapse  when  the  surrounding  fluid  is 
suddenly  exposed  to  a  pressure  increase.  He  derived  the 
continuity  and  momentum  equations  for  the  layer,  neglecting  the 
convective  and  dissipative  terms  and  assuming  that  the  volume 
fraction  of  gas  is  small  enough  to  warrant  such  approximations. 
However,  he  did  take  into  account  the  individual  bubble  radial 
motion  and  translation,  neglecting  viscous  effects.  Solving  the 
derived  system  of  equations.  Van  Wijngaarden  found  a 
considerable  increase  of  the  pressure  along  the  wall  due  to 
collective  effects. 

Recently  d'Agostino  and  Brennen  (1983),  using  an  approach 
similar  to  that  of  Van  Wijngaarden,  considered  the  unidimen¬ 
sional  problem  of  a  spherically  symmetric  cloud.  Considering 
small  harmonic  fluctuations  of  the  tar  field  pressure,  they 
investigated  the  cloud  response  for  various  oscillating 
frequencies.  Neglecting  bubble-liquid  relative  motion  and 
bubble  interactions,  they  obtained  very  interesting  results  on 
the  oscillations  of  the  whole  cloud  comprised  of  equal  sized 
bubbles.  Cases  for  which  the  natural  frequency  of  the  cloud  is 
much  smaller  than  the  bubble  frequency  were  found.  A  shielding 
of  the  cloud  core  by  the  outer  layer  was  obtained  when  the 
forcing  frequency  was  higher  than  the  individual  bubble 
frequency.  The  neglect  of  the  bubble  translation  velocity  and 
of  bubble  interactions  are,  however,  major  assumptions  which 
require  further  analysis. 

Mdrch  (1977,  1980),  concerned  with  ultrasonic  cavitation 
fields,  considered  the  collapse  of  a  hemispherical  "cluster" 
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near  a  wall,  which  he  extended  by  symmetry  to  the  case  of  a 
spherical  cloud.  He  characterized  the  cloud  only  by  its  radius 
and  a  uniform  void  fraction,  a,  constant  in  time,  and  developed 
the  following  model.  A  pressure  rise  in  the  liquid  leads  to  the 
formation  of  a  shock  wave  on  the  "cloud  boundary."  The  shock 
moves  toward  the  cloud  center  leaving  no  bubbles  behind  it  and 
thus  constitutes  the  cloud  boundary  at  each  time.  The  collapse 
time  of  a  spherical  cloud  was  found  to  be  a  times  the  Rayleigh 
collapse  time  of  a  spherical  bubble  of  the  same  initial  radius. 
Although  a  very  interesting  approach,  especially  for  the 
calculation  of  the  collapse  time,  this  model  (like  Rayleigh's 
model  for  spherical  bubble  collapse)  is  incapable  of  adequately 
calculating  the  pressure  field.  At  the  end  of  the  collapse  the 
cloud  radius  is  zero,  and  the  velocities  and  pressures  are 
infinite  since  the  model  does  not  allow  the  bubbles  to  contain 
noncondensibles.  The  main  physical  assumption  (presence  of  a 
shock  wave  dividing  the  space  in  two  regions,  one  containing 
bubbles  which  do  not  sense  the  pressure  variations  until  a  later 
time  stage,  and  another  one  where  all  bubbles  have  collapsed)  is 
valid  only  for  relatively  high  void  fractions.  The  case  of  a 
spherical  single  cavity  of  the  same  size  as  the  whole  cloud  is 
the  perfect  extreme  example  of  the  domain  of  validity  of  this 
approach.  Hanson  and  Mjbrch  (  1980)  and  M^rch  (1982)  extended  the 
same  model  to  a  cylindrical  cloud  and  a  layer  of  bubbles  on  a 
solid  plate. 

In  this  report  we  extend  the  singular  perturbation  theory 
earlier  developed  to  study  the  interaction  between  collapsing 
bubbles  (Chahine  and  Bovis,  1981)  to  an  asymptotic  theory  for 
the  collective  behavior  of  a  multibubble  system.  A  general 
theory  including  mass  and  heat  transfer  is  first  presented,  then 
within  the  constraints  of  simplifying  assumptions,  particular 
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studies  are  developed.  The  study  of  the  purely  dynamical 
behavior  of  the  bubble  cloud  is  the  most  advanced,  and  codes  for 
a  system  of  N  equi-sized  symmetrically  distributed  bubbles 
(Chahine  1982  a,b,c;  Chahine  and  Shen,  1985)  as  well  as  a  code 
for  any  given  distribution  of  bubbles  were  developed  and 
extensively  used.  For  instance,  the  symetrical  program  was 
applied  to  the  study  of  the  use  of  a  venturi  as  a  cavitation 
susceptibility  meter  and  has  been  recently  improved  to  include 
gravity  effects  and  the  corresponding  bubble  liquid  relative 
velocity.  Heat  (Chahine  and  Liu,  1983,  1984,  and  1985)  and  mass 
transfer  across  the  bubble-liquid  interface  are  the  subjects  of 
two  separate  studies.  In  both  cases  we  have  considered  the 
simplifying  assumptions  of  a  symmetrical  bubble  configuration 
and  of  a  thin  boundary  layer  in  which  either  temperature  or  gas 
concentration  varies.  Two  separate  codes  were  developed,  based 
on  these  two  studies,  and  some  illustrative  examples  are 
presented  here.  Liquid  compressibility  is  modeled  in  an 
additional  code  by  use  of  an  artificially  input  finite  sound 
speed  and  time  duration  for  propagation  from  one  bubble  to 
another.  Systematic  experimental  observations  of  spark 
generated  bubbles  near  complex  wall  geometries  simulating  a 
multibubble  system  were  made  in  order  to  validate  the  theory. 
Simultaneously,  pressuresgenerated  by  the  simulated  cloud  were 
measured  and  analyzed  (Chahine  and  Sirian,  1985).  Acoustic 
transmission  through  a  sheet  of  bubbles  was  investigated  using 
bubbles  generated  by  air  injection  and  emitting  and  receiving 
hydrophones . 


«  .  *  „  *  »  *  .  *  .  '  .  *  . 
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2.0  GENERAL  FORMULATION  OF  THE  PROBLEM 

As  a  first  step  to  study  the  general  problem  of  a  bubble 

cloud  in  a  flow  field  and  near  solid  boundaries,  let  us  consider 

a  cloud  of  bubbles  in  an  unbounded  medium  of  uniform  pressure, 

P*,  concentration  of  dissolved  gas,  C„,  and  temperature, 

T*,.  This  corresponds  to  the  case  where  the  size  of  the  cloud 

is  small  compared  to  the  flow  field  characteristic  length 

scale.  Poo,  Coo,  and  Too  are  then  the  local  values  of  the 

pressure,  the  concentration,  and  the  temperature  in  the  flow 

field  in  the  absence  of  the  cloud.  We  further  assume  the  liquid 

to  be  inviscid  and  incompressible  and  the  flow  irrotational . 

These  assumptions  are  commonly  accepted  and  are  justified  in 

cavitation  and  boiling  heat  transfer  studies  except  in  the  last 

phases  of  bubble  collapse.  The  neglect  of  finite  sound  speed 

effects  can  also  be  unacceptable  in  the  very  early  phases  of 

bubble  growth  (Baumeister  and  Hamill,  1969).  However,  we  will 

8 

not  be  concerned  with  these  early  times  (t<10  s) ,  especially 
since  numerical  experiments  have  shown  that  the  later  history  is 
influenced  very  little  by  the  value  of  the  initial  time  at  which 
the  computation  is  started.  The  bubble  cloud  behavior  is  sought 
when  the  ambient  pressure,  P«(t),  is  time  dependent. 

In  order  to  determine  the  flow  field  in  the  bubble  liquid 
medium  and  to  obtain  the  motion  and  deformation  of  any  bubble  in 
the  cloud,  one  must  solve  the  Laplace  equation  for  the  velocity 
potential,  * 

At  =  0  ,  (2.1  ) 


subjected  to  kinematic  and  dynamic  conditions  on  the  bubbles' 
surfaces . 
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The  equation  of  a  bubble  surface  in  a  coordinate  system 
moving  with  velocity  6*-  in  the  direction  ez ,  is 
r  *  Ri(9,<fr,t).  and  n*  are  respectively  the  local 

curvature  of  the  surface  of  bubble  B1  and  its  unit  normal 
vector  at  the  point  M(r,9,$).  y  and  Tr  are  respectively  the 
surface  tension  of  the  liquid  and  its  temperature  at  the  bubble 
wall.  Pg  and  Pv  are  the  partial  pressures  of  the 
noncondensible  gas  and  the  vapor  inside  the  bubble.  The 
boundary  conditions  can  then  be  written: 

V*  •  n1!  =  [R1  e  +  6le  ]  n1  ,  (2.2) 

r=R1( 9 , $,t) 


p  [  $  -  bxe  +  i  I  V*|  2  ] 

2  1  1  r=R1  {  9 ,  $,  t ] 


=  P  (t) 


-  Pt<V  ' 


p‘(ti 


+  2  yMTrJ^M  9,  *,t)  ,  (2.3) 


where  $  and  the  operator  7  are  expressed  in  the  moving 
coordinate  system.  Due  to  the  low  value  of  the  vapor  density, 
pv,  the  pressure  of  the  gas  and  the  vapor  inside  the  bubble 
can  be  assumed  to  be  uniform  as  long  as  the  spherical  symmetry 
is  preserved. 


Due  to  the  relatively  very  short  time  of  vaporization 
compared  to  bubble  dynamics  and  gas  diffusion  characteristic 
times,  the  vapor  is  considered  to  instantaneously  flow  in  and 
out  of  the  bubble,  and  Py  is  taken  equal  to  the  value  of  the 
equilibrium  vapor  pressure  of  the  liquid  at  the  bubble  wall 
temperature.  When  the  bubble  shape  deviates  moderately  from  a 
sphere  we  will  assume  that  both  the  temperature  along  the  bubble 
wall  and  the  value  of  the  vapor  pressure  vary  accordingly. 

Under  this  assumption,  the  pressure,  Pv,  may  be  uniform  inside 
the  bubble  far  from  the  bubble  surface  but  accomodates  itself  to 
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the  temperature  controlled  value  in  the  vicinity  of  the 
interface.  More  details  of  the  way  this  happens  inside  the 
bubble  are  not  needed  here  since  the  flow  field  of  the  vapor  is 
of  no  relevance  as  long  as  the  velocities  are  subsonic.  The  gas 
and  the  vapor  are  taken  to  comprise  an  ideal  gas  mixture  and 
possess  the  same  temperature.  The  partial  pressure  of  the  gas, 
Pg,  is  determined  using  the  ideal  gas  law  as  well  as  an  energy 
balance  of  the  bubble  content.  The  concentration  of  dissolved 
gas  at  the  bubble  surface  is  related  to  the  partial  pressure  of 
the  gas  within  the  bubble  by  Henry's  Law, 

Pg  =  H  C 1 ( R1 )  ,  (2.4) 

where  H  is  Henry's  constant. 

The  values  of  the  equilibrium  vapor  pressure  PV(TR)  and 
of  the  surface  tension  y(Tr)  constitute  the  coupling  between 
the  dynamic  and  the  heat  problems. 

Similarily  the  value  of  Pg(t)  constitutes  the  coupling 
between  the  dynamic  and  the  mass  transfer  (gas  diffusion) 
problem.  To  determine  the  temperature  at  the  bubble  wall, 

TR( 0 , $, t ) ,  one  needs  to  solve  the  energy  equation  in  the 
liquid , 

T  +  7*  •  VT  =  D  •  AT,  (2.5) 

where  D  is  the  thermal  diffusivity  of  the  liquid.  Equation 
(2.5)  is  subjected  to  a  boundary  condition  on  the  bubble  wall 
stating  that  the  heat  locally  lost  at  any  point  of  the  interface 
is  used  to  vaporize  an  amount  of  liquid  determined  by  the  local 
bubble  volume  expansion  rate.  If  pv  is  the  vapor  density,  L 


Tracor  Hydronautics 


-9- 


the  latent  heat  of  the  liquid,  and  K  its  thermal  conductivity, 
the  heat  balance  equation  over  the  bubble  surface  can  be  written 
in  spherical  coordinates: 


K  II 

dn1 


r=R  ( 0, $,t) 


i  ^ 

R  {  0,  4>,t)  sin  0  d  9d  <#»  = 


=  PvL  ^  [//  R  sinfl  de  d*] 


(2.6) 


This  equation  is  satisfied  if  the  following  elementary 
equilibrium  equation  applies  locally  at  the  bubble  surface: 


9T 

9nj 


pvL 


R 


r=R  (9, *,t) 


(2.7) 


The  concentration  of  gas  in  the  liquid,  C(r,0,$,t),  is 
determined  by  the  diffusion  equation 

C  +  •  _VC  =  Dg  •  AC  ,  (2.8) 


where  Dg  is  the  molar  diffusivity  of  the  gas  in  the  liquid. 
This  equation  is  subjected  at  the  bubble  wall  to  the  condition 
that  the  concentration  at  the  interface  is  equal  to  the 
saturation  concentration  at  the  partial  pressure  of  the  gas  in 
the  bubble 

C(r1)  =  C„f  .  (2.9) 


Csat  is  related  to  Pg  by  Henry's  law,  (2.4).  The  molar  rate 
of  transport  of  gas  across  the  interface,  ng ,  is  related  to 
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the  gradient  of  gas  concentration  in  the  liquid  at  the  bubble 
wall  by  the  relation 

j.2 

R  (9,$,t)  sm  0  d0  d$  (2.10) 

r=R1  (  0 ,  t ) 

Time  integration  of  (2.10)  determines  at  every  instant  the  total 
number  of  moles  of  gas,  ng ,  in  the  bubble.  The  partial  gas 
pressure  inside  the  bubble  is  then  related  to  its  volume  and  to 
ng  by  a  perfect  gas  equation  of  state.  The  only  remaining 
variable  needed  to  apply  this  perfect  gas  law  is  the  gas 
temperature  inside  the  bubble.  The  application  of  the  first  law 
energy  balance  for  the  bubble  interior  relates  this  temperature 
to  the  boundary  work  due  to  volume  change,  the  enthalpy  transfer 
due  to  transport  of  gas  and  vapor  across  the  bubble  wall,  and 
change  in  internal  energy  of  the  bubble  contents.  This  is 
expressed  in  detail  in  Section  6. 

Equations  (2.1)  to  (2.10)  form,  with  the  bubble  energy 
equation  and  the  initial  and  at-infinity  conditions  (known  T*, 
and  P.(t)),  a  complete  set  of  equations  which  must  be  solved 
to  determine  the  flow,  concentration,  and  temperature  fields. 


n  =  D 

g  g 


// 


3C 

3n] 
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3.0  DYNAMICS  PROBLEM:  SINGULAR-PERTURBATION  THEORY 


3.1  Method 


The  following  approach  is  applicable  to  a  cloud  of  bubbles 
of  low  void  fraction.  Provided  that  the  characteristic  size  of 
a  bubble  in  the  cloud,  rb0»  is  small  compared  to  the 
characteristic  distance  between  bubbles,  Iq,  we  can  assume,  in 
the  absence  of  an  initial  relative  velocity  between  the  bubble 
and  the  surrounding  fluid,  that  each  of  the  individual  bubbles 
reacts,  in  first  approximation,  to  the  local  pressure  variations 
spherically  as  if  isolated.  To  the  following  order  of 
approximation,  interactions  between  bubbles  induce  bubble  motion 
and  deformation  and  are  taken  into  account.  This  approach  is  an 
extension  of  the  earlier  studies  by  Chahine  and  Bovis  (1981)  and 
Bovis  and  Chahine  (1981)  on  the  collapse  of  a  bubble  near  a 
solid  wall  and  a  free  surface  (later  presented  more  generally 
for  nonspherical  bubbles  by  Chahine,  1982). 

Since  the  problem  possesses  two  different  length  scales, 
lQ  and  rbQ»  we  can  consider  two  subproblems:  one 
concerned  with  the  macroscale  and  the  other  one  with  the 
microscale.  The  "outer  problem"  is  that  considered  when  the 
reference  length  is  set  to  be  iQ.  This  problem  addresses  the 
macrobehavior  of  the  cloud,  and  the  bubbles  appear  in  it  only  as 
singularities.  The  "inner  problem"  is  that  considered  when  the 
lengths  are  normalized  by  rb0  and  its  solution  applies  to 
the  vicinity  of  the  considered  individual  bubble  of  center 
B*.  The  presence  of  the  other  bubbles,  all  located  at 
infinity  in  the  "inner  problem",  is  sensed  only  by  means  of  the 
matching  condition  with  the  "outer  problem".  That  is  to  say, 
physically  the  boundary  conditions  at  infinity  for  the  "inner 
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problem"  are  obtained,  at  each  order  of  approximation,  by  the 
asymptotic  behavior  of  the  outer  solution  in  the  vicinity  of 
B*-.  Mathematically,  one  has  to  match  term  by  term  the  inner 
expansion  of  the  outer  solution  with  the  outer  expansion  of  the 
inner  solution  using  the  same  asymptotic  sequence  in  the  two 
expansions . 

3 . 2  First  Order  Approximation 

The  determination  of  the  flow  field  and  the  dynamics  of  any 
of  the  individual  bubbles,  B^,  is  accessible  once  the  boundary 
conditions  at  infinity  in  the  corresponding  "inner  region"  are 
known.  Here  we  impose  the  restrictive  assumption  that  the  void 
fraction  is  low  enough  so  that  the  information  about  the 
variation  of  the  ambient  pressure  around  the  cloud,  P„(t),  is 
transmitted  to  the  microscale  in  a  time  scale  much  shorter  than 
the  bubble  collapse  time.  Therefore,  in  the  absence  of  a  slip 
velocity  between  the  considered  bubble  and  the  surrounding  fluid 
and  when  interactions  are  neglected,  the  only  boundary  condition 
at  infinity  is  the  imposed  pressure  variation  P„(t). 

When  a  finite  value  of  the  sound  speed  in  the  cloud  medium 
must  be  accounted  for,  P*(t)  can  be  approximated  by  Pffl(t' ) 
where  t'  *  t  -x/c  and  x  represents  the  minimum  thickness  of  the 
cloud  shell  at  the  bubble  location.  To  address  the  problem  more 
rigourously  one  must  solve  the  two-phase  medium  flow  equation 
using,  for  example,  a  method  similar  to  that  proposed  in  Section 
7.0.  The  "inner  problem"  is  therefore  spherically  symmetric  and 
its  solution  is  given  by  the  well  known  Rayleigh-Plesset 
equation.  With  the  assumption  that  the  liquid  is  inviscid  and 
incompressible,  this  equation  can  be  written  as  follows  when 
mass  and  heat  transfer  are  neglected 
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a  ^ 

o  o 


3 

1 


^  =  -P 
o  1 


(t)  +  P„  (a 


-  3k 


-  1)  +  W-1  (1  - 


o 


(3.1  ) 


This  equation  is  slightly  modified  when  P7  and  Pg  are 
dependent  on  the  transfer  of  heat  or  mass  across  the  surface. 
It  is  rewritten  in  those  cases  in  the  corresponding  sections. 
In  this  equation,  where  the  superscript  i  is  omitted  for 
convenience,  a0(t)  is  the  radius  of  bubble  B1  normalized  by 
rt>0*  The  times  are  normalized  by  the  Rayleigh  time,  T0, 
based  on  rbQ  and  a  characteristic  value  of  the  pressure 
variations,  AP.  All  pressures  are  normalized  by  AP.  PgQ  is 
the  initial  normalized  gas  pressure  in  the  bubble.  Thus, 


T 


o 


p'/2/< 


AP) 


1/2 


(3.2) 


P  =  P  /(AP)  ,  (3.3) 

yo  yo 

P  Jt)  =  (PJt)  -  PQ)  /  (AP)  .  (3.4) 


The  Weber  number  is  related  to  the  surface  tension, 

rb0  by; 


We  =  rb  •  (  AP)  /  2y  . 
o 


y,  AP,  and 


(3.5) 


The  noncondensible  gas  pressure  inside  the  bubble,  Pa ,  is 

3k  y 

assumed  to  have  a  polytropic  behavior,  Pg  ao  =  constant, 
where  k  is  the  polytropic  coefficient  (1  <  k  <  Cp/cv). 

For  a  given  P»(t),  Equation  (3.1)  can  be  solved  for  the 
variation  of  the  bubble  radius,  aQ(t).  This  allows  the 
subsequent  determination  of  the  pressure  field  around  the  bubble 
of  center  B1,  by  the  use  of 
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V8*' 


t)  =  P  (t)  +  r 


-  l 


(2  a. 


+  a. 


ao)  " 


a*/2r\(3.6) 


where  r  is  the  distance  between  B1  and  a  given  point  M  in  the 
fluid . 


3.3  Interactions 


When  interactions  cannot  be  neglected,  but  an  "inner 
region"  enclosing  the  bubble  B1  can  still  be  defined,  the 
boundary  conditions  at  infinity  can  be  much  more  complex. 

First,  the  macroscale  pressure  in  the  cloud  at  B*-,  P(B^,  t)  , 
can  be  very  different  from  the  imposed  far  field  pressure 
Poolt).  Second,  a  relative  velocity  between  the  bubble  and  the 
surrounding  fluid,  U(Bi,  r,  t)  can  exist  causing  the  bubble  to 
be  nonspherical .  Both  P  and  U  can  be  determined  only  by  solving 
the  equations  of  motion  of  the  two-phase  medium  (see  Section 
7.0).  In  this  section  we  will  limit  ourselves  to  a  small 
perturbation  theory  whose  interest  will  be  to  give  the  behavior 
of  the  solution  when  the  perturbation  grows  continuously.  In 
that  case  PfB*-,  t) ,  which  is  the  driving  pressure  for  the 
collapse  of  the  bubble  B1,  is  only  a  perturbation  of  the 
imposed  far  field  pressure,  P»(t) ,  and  U(B^,  r,  t)  is  a 
perturbation  of  the  spherical  velocity  due  to  the  bubble  volume 
variation. 

If  we  assume  that  the  liquid  flow  is  irrotational ,  we  can 
define  a  velocity  potential  for  the  macroscale  ("outer 
problem"),  ${BX,  t) ,  and  a  velocity  potential  for  the 
microscale  ("inner  problem"),  *i(B1M,  t)  such  as: 
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The  matching  condition  between  these  two  potentials  expresses 
the  at-infinity  conditions  for  a1,  and  replaces  the  conditions 
on  P  ( B* ,  t)  and  U(B^,  r,  t) .  This  is  written: 

*  .  . 

lim^  <|>1(B1M,  t)  =  0  +  e  ♦1(B1,  t)  +  e2*2(B  ,  t)  +  ... 

BXM  *  »  (3.8) 


Since  bubble  interactions  vanish  when  e  goes  to  zero,  in  the 
absence  of  relative  velocity  between  the  surrounding  fluid  and 
the  bubbles  in  the  cloud,  $0  =  0.  *lf  *2'  •••  are  the 
contributions  of  the  whole  cloud  to  the  boundary  condition  at 
infinity  for  the  inner  problem  (i).  Using  the  results  obtained 
with  the  interaction  of  two  bubbles  (Chahine  and  Bovis,  1981) 
and  the  property  of  addition  of  potential  flows,  this  condition 
can  be  written: 


where  the  superscript  (j)  denotes  quantities  corresponding  to 
the  other  bubbles,  B3 .  is  the  initial  distance  between 

the  bubble  centers  B*-  and  B3  .  9^3  is  the  angle  MB^B3  , 

and  r  the  distance  B*M,  where  M  is  a  field  point  in  the  fluid 
(see  Figure  3-1  )  .  (cose)  is  the  Legendre  polynomial  of 

order  n  and  argument  cose,  q^  is  the  correction  of  order  en 

of  the  strength,  q3  =  a3  (a3)2,  of  the  source  representing 
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the  first-approximation  spherical  oscillations  of  the  bubble 
B3  . 


Expressed  in  physical  terms  (velocities,  pressures),  the 

boundary  condition  (3.9)  states  that  the  first  order  correction, 

(e),  to  the  nonperturbed  spherical  behavior,  aQ(t),  of  the 

bubble  B*  is  a  spherical  modification  of  the  collapse  driving 

pressure.  This  would  introduce,  as  for  two  bubbles,  a  spherical 

correction  ai(t)  of  the  variations  aQ(t).  At  the  following 

order,  (e2),  a  second  correction  of  the  uniform  pressure 

appears,  as  well  as  a  uniform  velocity  field  accounting  for  a 

slip  velocity  between  the  bubble  and  the  surrounding  fluid. 

Again,  as  in  the  two-bubble  case,  this  induces  a  spherical 

correction,  a2(t),  of  a^t),  and  a  nonspherical  correction 

f^(t)  cos  0*9 ,  where  0*9  is  a  direction  to  be  determined  from 

all  the  9^3.  Things  become  more  complex  at  the  order  of 

expansion  e3,  where  in  addition  to  the  uniform  pressure  and 

i  i  ■ 

velocity  corrections,  a3(t)  and  f3(t)  cos  0*9 ,  a  velocity 
gradient  is  to  be  accounted  for,  to  generate  a  nonspherical 
correction,  g3(t)  (3  cos 20*9-l )/2 . 

As  a  result,  the  equation  of  the  surface  of  the  bubble  B1 
can  be  written  in  the  form: 

R  01*,  t)  =  a*(t)  +  e  a  i(  t)  +  e2  [a2(t)  +  f2(t)  •  cos01;^|  + 

+  e3  [a  3(  t)  +  f3(t)  cos0lj  +  g3(t)  (cos0lj)]  +  ...,  (3.10) 

where  a^t)  is  given  by  the  Rayleigh-Plesset  Equation  (3.1), 
while  the  other  corrections  are  obtained  by  solving  the 
following  differential  equations  in  which  the  superscripts  (i) 
are  omitted  for  convenience: 
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a0a  i 


3aoa  i 


+  aiFo(ao' 


f 


a  a  2  +  3a  i  2  + 
0  *  0  £ 

a2F0  +  F  l(a0'  a  1'  We' 

P  ,  K) 
^0 

a  a  3  +  3a  a  3  + 
0  3  0  3 

a3F0  +  F2(aQ»al»a2  'We 

,P  ,K> 

’O 

-•>-(?)  «• 

••  .  • 

a  d?  +  3a  do  = 
0  *■  0  *■ 

+  a  <q3) 
00 

,  (3.11) 

••  .  . 

a  d,  +  3a  d3  + 
0  3  0  3 

3F  3  ( aQ  ,  a  ! )  d  2  =  - 

W3> 

2 

1  (®oq'+ao’>+Pqi)  ' 

ao^3  +  3aog3  - 

(ao_6/Weao)  g3  =  Ej  " 

V‘o  > 

1  (ao£io+2aoa'o<5o)  • 

In  these  equations  F0,  Flf  F2,  F3  are  known  functions 
depending  on  the  physical  constants,  We  and  Pg0,  and  on 
the  calculated  preceding  orders  of  approximation.  The 
deformations  f 3  of  the  bubble  B1  and  the  motion  of  its 

center  toward  B3 ,  l2,  £3,  have  been  replaced  by  d2,  d3  which 
indicate  the  total  motion  of.  the  point  (Figure  3.1). 

d  2  =  f  2  ~  l2i  ^  3  =  £3  ~  *3  (3.12) 

When  all  the  initial  radii  of  the  bubbles  in  the  cloud  are 
identical,  these  right-hand  sides  are  obtained  by  multiplying 
the  two-bubble  case  right-hand  sides  by  one  of  the  geometrical 
constants  Cj,  c2,  c3: 
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ci  =  Ej  (  lo/lo3)  ' 

c2  =  (  V4j)2  *  cos0lj  ,  (3.13) 

C  3  =  Ej  (  £  j  )  3  *  <^2  (COS0lj) 

We  can  now  compute  the  behavior  of  B1  by  solving  the  obtained 
differential  equation  using  a  multi-Runge-Kutta  procedure.  The 
behavior  of  the  whole  cloud  can  then  be  obtained. 


Comparison  of  Equations  (3.11)  with  those  obtained  in  the 
case  of  two-bubbles  shows  that  the  N-bubbles  in  the  cloud  other 
than  can  be  replaced  by  a  unique  bubble  of  strength,  q|9 , 
located  at  G^,  a  distance  1^3  from  B^  in  the  direction 
defined  by  the  angle  MB^-G*  =  0*D  .  As  this  equivalent 
bubble  should  induce  the  same  pressure  and  velocities  as  defined 
by  (3.9),  its  characteristics  are  obtained  by  the  equations: 
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where  e^g  and  e^j  are  unit  vectors  of  the  directions 
and  B^bD  respectively  (Figure  3.1),  and  n  is  the 
order  of  approximation.  These  equations  define  the  angle  0^9, 
and  the  direction  in  which  dn(t)  is  measured  (Equations  (3.10) 
and  (3.12)). 
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3.4  Pressure  Field 


Once  the  "inner  problems"  are  solved,  the  nondimensional 
outer  potential,  *,  can  be  written: 


♦(M,t) 


cos  0*^+  0  (  e  3)  , 
(3.16)-* 


where  bars  denote  nondimensional  "outer"  quantities,  and  tildes 
nondimensional  "inner"  quantities. 
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T0  is  the  characteristic  time  of  the  bubble  collapse  and  r^ 
is  the  distance  between  a  field  point  M  and  B*.  The  Bernoulli 
equation  enables  one  to  calculate  the  pressure,  P,  using 
Equation  (3.16).  In  the  nondimensional  form  we  have: 


p(M,t) 


p(M,t)  -  P«(t) 
AP 


ii  1 

at  2 


(3.18) 


AP  is  the  amplitude  of  the  pressure  driving  the  collapse  and  t  = 
t/T0 ,  where 

To  =  rb  V p/  AP '  *  (3.19) 

o 


As  an  illustration,  in  a  uniform  field  of  bubbles  any  bubble  has 
the  same  geometrical  position  relative  to  the  others,  and  thus 
the  same  behavior.  The  general  expression  (3.8)  simplifies 
considerably  to  become: 
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p(M,t)  =  (e  qQ  +  £2  q!  +  e3  q2  +  e4  q3)  Ei  + 
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In  this  expression,  the  two  first  summations  are  geometrical 
constants  similiar  to  c3,  c2  (3.4).  The  last  one  is  more 
complex,  but  is  more  easily  calculated  when  written  as  follows 

-i  2  —  -  2 


’!ti  ?]  ■  h  (ih)  e*im] ' 


(3.21  ) 


where  eim  is  the  unit  vector  of  the  direction  BXM.  If  one 
knows  the  direction,  MV0 ,  of  the  velocity  at  M,  at  the  first 
order  of  approximation,  and  if  axv  is  the  angle  BXMV0 
( Figure  3.1),  then : 


/  "  (?)] 2 


(3.22) 


In  the  general  case  when  the  bubbles  do  not  all  behave 
identically  an  expression  similar  to  (3.20)  applies  in  which  the 
summation  sign  concerns  not  only  the  geometry  but  also  the 
source  and  dipole  intensities,  q1  and  h1. 


Tracor  Hydronautics 


-21  - 


4.0  DYNAMICS  PROBLEM:  NUMERICAL  STUDY  AND  RESULTS 
4 . 1  Numerical  Codes  Developed 

Several  numerical  codes  have  been  developed  during  this 
study.  For  treating  only  dynamics  (mass  and  heat  transfer 
neglected),  two  main  numerical  programs  can  be  distinguished: 
CLDSYM  and  MULTIBBL.  Other  programs  were  also  implemented  that 
are  improvements  or  modifications  of  these  two  codes.  The  two 
main  programs  are  distinguished  by  the  cloud  geometry:  either 
symmetrical  or  arbitrarily  described. 

When  all  the  bubbles  have  the  same  initial  size,  and  when 
their  initial  geometrical  configuration  is  symmetrical,  they  all 
behave  identically  by  reason  of  symmetry.  At  a  given  order,  all 

the  functions  q3 ,  which  appear  in  Equation  (3.9)  are  the 

n 

same.  Therefore,  the  summation  applies  only  to  geometrical 
constants  which  are  known  at  the  beginning  of  the  computation 
and  depend  only  on  the  bubbles'  distribution.  The  three 
constants,  clr  c Zt  and  c3,  defined  in  (3.13)  are  then  needed  to 
compute  the  seven  components  of  R(  0,t)  up  to  the  order  e3 
included.  They  are  used  as  input  to  the  numerical  program 
dubbed  CLDSYM,  which  uses  a  multi-Runge-Kutta  procedure  to 
compute  the  bubble  motion  as  well  as  the  generated  flow  and 
pressure  field.  Figure  4.1  shows  a  simplified  flow  chart  of  the 
numerical  procedure. 

When  the  bubbles  have  different  initial  radii,  and  when 
their  distribution  in  the  liquid  does  not  correspond  to  any 
particular  symmetrical  configuration,  the  simplification 
described  above  does  not  apply.  A  second  numerical  program 
dubbed  MULTIBBL  has  been  developed  for  this  more  general  case. 
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The  flow  chart  in  Figure  4.2  describes  this  program.  In  this 
flow  chart  Qm  represents  the  right-hand  side  of  Equation  (2), 
while  qm  represents  the  strength  of  the  flow  sources  obtained 
at  the  order  of  approximation  m.  This  program  allows  the 
computation  of  more  realistic  cases  where  the  bubble  sizes  and 
distribution  are  input  and  do  not  correspond  to  any  idealized 
perfect  configuration  as  with  CLDSYM.  More  realistic  pressures 
in  the  far  field  as  well  as  on  nearby  walls  are  then  obtained 
since  the  different  bubble  collapses  are  not  in  phase  as  for  the 
case  of  a  symmetrical  configuration  of  equal  sized  bubbles. 

A  third  program  dubbed  CLDMAIN2  was  implemented  for  the 
cases  where  the  computation  spans  over  several  thousand  times 
the  bubble  characteristic  time.  This  program  adapts  a  variable 
time  step  scheme  to  the  symmetrical  bubble  cloud  case.  An 
England  Runge-Kutta  procedure  is  used.  This  method  computes  the 
numerical  errors  introduced  at  each  time  step  and  determines  the 
time  increment  needed  to  both  minimize  computation  time  and  to 
not  exceed  an  imposed  acceptable  margin  of  error.  A  major 
advantage  of  this  code  is  its  capacity  to  describe  very  rapid 
bubble  radius  variations  such  as  occurs  during  strong  explosive 
growth  or  collapse.  This  code  was  used  for  the  analysis  of 
bubble  dynamics  in  a  venturi  cavitation  susceptibility  meter  and 
about  a  hydrofoil. 

A  fourth  program  dubbed  MULTICOMP  was  derived  from  MULTIBBL 
and  employs  a  finite  sound  speed  input.  The  information 
propagation  from  one  bubble  to  another  is  retarded  by  an  amount 
equal  to  the  ratio  between  the  interdistance  and  this  finite 
sound  speed.  This  program  was  used  to  study  the  influence  of 
time  lags  creating  phase  differences  between  bubble  oscilla¬ 
tions  . 


Tracor  Hydronautics 


-23- 


These  four  programs  were  used  to  analyze  several  variables. 
The  objectives  of  exercising  these  programs  were  to  debug  them, 
demonstrate  their  performance,  and  obtain  trends  of  behavior. 
Below  we  illustrate  some  of  the  results  obtained  in  order  to 
gain  some  understanding  of  the  behavior  of  a  bubble  cloud  and 
how  it  differs  from  the  behavior  of  isolated  bubbles.  The 
influence  of  the  main  parameters  will  be  stressed. 

4 . 2  Preliminary  Discussion  on  Bubble  Interaction 

It  is  helpful,  in  order  to  understand  the  various  numerical 
results,  to  consider  the  pressure  field  generated  by  the 
collapse  of  an  individual  bubble.  Let  us  consider  the  pressure 
history  at  a  point  M  at  a  distance  Z0  from  an  isolated  bubble 
when  the  ambient  pressure  is  submitted  to  a  sudden  pressure 
jump.  As  we  can  see  from  Figure  4.3,  the  resulting  perturbation 
pressure,  i.e.,  the  difference  between  the  pressure  at  £0  and 
the  far  field  pressure,  is  negative  for  t<0.75  for  the 
considered  case.  Later,  for  t>0.75,  the  pressure  at  £0 
increases  rapidly  to  2.3  times  the  imposed  pressure  at  the  end 
of  the  bubble  collapse.  This  observation  is  a  key  to  the 
understanding  of  multibubble  behavior  in  a  sudden  pressure 
increase.  Indeed,  if  a  fictitious  bubble  which  does  not 
influence  the  process  is  centered  in  B  lf  it  will  first  sense  a 
less  important  and  more  gradual  increase  in  the  surrounding 
pressure  at  the  beginning  of  the  collapse.  In  the  case 
considered  (Figure  4.3),  instead  of  a  nondimensional  surge  of 
the  pressure  from  0  to  1 ,  P  jumps  only  to  0.84  then  rises 
slowly,  not  attaining  1  until  t=0.75.  This  induces  a  relative 
decrease  in  the  intensity  of  the  first  phase  of  the  bubble 
collapse  and  affects  the  bubble  radius  variation  with  time 
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exactly  as  seen  later  in  Figure  4-6.  Later,  the  sharp  rise  in 
the  pressure  field  makes  the  subsequent  end  of  the  collapse  much 
more  violent. 

Figure  4-4  illustrates  the  case  where  the  ambient  pressure 
field  is  modified  by  the  collapse  of  more  than  one  bubble.  This 
figure  is  intended  to  explain  both  cumulative  effects  due  to  a 
multibubble  system  and  modification  due  to  different  bubble 
size,  bubble  interdistance,  and  delay  times  (compressibility) 
effects.  The  figure  illustrates  these  effects  by  considering 
two  bubbles.  Had  the  two  bubbles,  A  and  B,  been  of  the  same 
size  and  at  equal  distance  from  M,  the  two  corresponding 
pressure  time  profiles  would  have  superimposed. 

In  this  case,  both  effects  described  above  (initial  relative 
pressure  reduction  and  subsequent  dramatic  pressure  increase) 
will  be  doubled  (without  any  additional  consideration  of  the 
interaction  such  as  modification  of  the  pressure  histories) . 

The  interaction  between  bubbles  A  and  B  would  amplify  the 
phenomena  described  above.  Now,  if  bubbles  A  and  B  are  either 
at  a  different  distance  from  M  or  are  not  of  equal  size,  the 
pressure  profiles  will  not  overlap.  If  they  are  of  equal  size 
but  at  different  distances  from  M,  the  pressure  peaks  would 
occur  simultaneously  unless  compressibility  effects  are 
included.  A  delay  time  corresponding  to  the  difference  in  the 
travel  paths  will  then  occur  between  the  two  peaks. 

From  the  above  discussion  it  is  evident  that  the  behavior 
of  a  fictitious  bubble  located  at  M  will  strongly  depend  on  both 
the  sizes  of  the  other  interacting  bubbles  and  their  space 
distribution.  Compressibility  would  have  the  same  effect  as  a 
nonsymmetric  bubble  distribution  in  addition  to  its  influence  on 
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the  collapse  intensity  and  amplitude  of  the  generated  pressures. 
It  would  therefore  tend  to  reduce  the  collective  effect  expected 
from  the  incompressible  symmetrical  case. 

4 . 3  Simple  Case  of  a  Step  Function  for  the  Ambient  Pressure 
Field 


4.' 


We  consider  a  distribution  of  bubbles  centered  on  the 
surface  of  a  sphere,  and  we  admit  that  each  of  the  bubbles  has 
the  same  position  relative  to  the  others.  In  this  case  the 
numerical  computation  time  is  reduced.  As  examples,  we  will 
consider  the  bubble  behavior  and  the  pressures  generated  for  two 
types  of  ambient  pressure  time  functions.  In  both  cases  the 
bubbles  are  at  equilibrium  with  the  ambient  pressure,  PQ,  at 
t  =  0.  Then,  in  case  A,  the  ambient  pressure  jumps  to  a  new 
constant  value,  PQ  +  AP,  at  the  following  instant.  In  the 
second  case  B,  the  imposed  pressure  drops  first  to  a  constant 
value,  PQ  -  AP,  keeps  this  value  until  t  =  AT,  and  then  comes 
up  again  to  the  initial  value  P0  (Figure  4.5).  As  an 
illustration  we  will  consider  the  pressures  generated  a)  in  the 
center  of  the  sphere;  b)  at  the  location  of  the  bubble  B1  if 
it  were  removed  and  c)  at  a  point  outside  the  cloud  at  a 
distance  rbQ  from  B1.  We  will  compare  the  results  with 
the  isolated  bubble  case. 


Knowing  the  initial  bubble  configuration  and  thus  c)f  c2, 
and  C3  the  relation  between  the  cloud  radius,  R,  and  1q  is; 

R  =  1/2  £0  c1/c2.  In  the  cloud  center,  position  (a),  the 
three  summations  in  (3.20)  have,  respectively,  the  values  N/R, 
N/R,  and  0.  In  position  (b)  these  values  are  clf  c2,  and  c3, 
and  at  a  distance  rt>0  from  B1  the  values  are  approximated 
by  C[  +  e"1,  c2  -  e~ 2 ,  and  (c2  +  e-2)2. 


Tracor  Hydronautics 


-26- 


Various  spherically  symmetrical  cloud  configurations  were 
investigated  numerically.  In  Figure  4.6  the  results  of  five 
different  computations  for  a  sudden  jump  in  the  imposed  ambient 
pressure  are  compared,  expansions  being  conducted  up  to  e3.  The 
ratio,  e  =  r5o/£0,  was  kept  constant  and  at  a  value  of 
0.05.  The  cases  of  two,  three,  and  twelve  bubbles  are  presented 
together  with  that  of  an  isolated  bubble.  The  fifth  case  is  an 
intermediate  situation  between  the  configurations  of  three  and 
twelve  bubbles.  This  case  is  arbitrary  and  is  only  determined 
by  the  choice  of  ci,  c2,  and  C3.  In  each  case  the  variation 
with  time  of  the  distance,  B1E1,  between  the  extreme  point 
on  a  bubble  E*,  and  its  initial  center,  B* ,  is  chosen  to 
represent  the  bubble  dynamics.  Taking  the  bubble  collapse  in  an 
unbounded  fluid  as  reference,  it  is  easy  to  see  from  Figure  4.6 
how  increasing  the  number  of  bubbles  changes  the  dynamics  of  the 
one  studied.  We  can  observe  first  that,  during  the  early  slow 
phase  of  the  implosion  process,  the  collapse  is  significantly 
delayed.  At  any  given  nondimensionsal  time  the  distance  between 
B*  and  E*  (and  simultaneously  the  bubble  characteristic 
size)  is  greater  when  the  number,  N,  of  interacting  bubbles 
increases.  Then,  in  the  final  phase  of  the  implosion  the 
tendency  is  reversed:  the  phenomenon  speeds  up  and,  in  a 
shorter  total  implosion  time,  the  final  velocities  of  the  motion 
are  higher  when  N  increases.  The  effect  is  explained  by 
accounting  for  the  modification  of  the  driving  pressure  of  the 
collapse  of  any  bubble  due  to  the  dynamics  of  the  other  bubbles 
as  described  in  Section  4.2. 

Figure  4.7  shows  the  behavior  of  the  bubbles  in  the  case  of 
a  pressure  variation  of  type  B  (Figure  4.5).  The  cases  of  an 
isolated  bubble  and  two,  three,  five,  and  twelve  bubbles  are 
investigated  again,  and  the  variations  of  B^E1  with  time  are 
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represented.  The  ratio  e  and  the  duration  AT  of  the  pressure 
drop  are  kept  constant  and  at  the  particular  values  of  0.1  and 
0.8  respectively.  Here,  as  in  the  preceding  figure,  noticeable 
changes  can  be  observed  when  the  degree  of  interaction 
increases.  First,  the  growth  is  slowed  down  and  retarded  in 
comparison  with  the  isolated  case.  Then,  the  collapse  is 
accelerated  and  as  a  result  the  total  implosion  time  decreases 
with  an  increase  in  the  number  of  bubbles,  N.  While  for  N  =  2, 
the  total  implosion  time  is  greater  than  that  of  an  isolated 
bubble,  for  N  =  12  the  time  is  significantly  smaller.  As  we 
will  see  below  this  acceleration  of  the  collapse  makes  the 
generated  pressures  at  the  end  of  the  collapse  higher  than  for 
the  single  bubble  case. 

Figure  4.8  compares  for  the  same  cloud  configuration 
(twelve  bubble,  e  =  0.1)  the  bubble  behavior  for  three  values  of 
the  duration,  AT,  of  the  pressure  drop.  The  greater  AT  is,  the 
longer  the  bubble  is  allowed  to  grow.  As  a  result  the  maximum 
size  it  attains  is  bigger,  but  its  lifetime  is  smaller.  Thus, 
the  resulting  collapse  is  much  stronger. 

As  explained  in  Section  4.2,  using  Figure  4.3,  the 
observations  made  above  for  a  sudden  pressure  jump  can  be 
explained  by  considering  how  the  presence  of  a  reacting  bubble 
modifies  the  local  pressure  field  from  that  imposed  at 
infinity.  The  same  type  of  reasoning  can  be  applied  to  the  case 
of  a  finite-time  pressure  drop.  As  we  can  see  in  Figure  4.9,  in 
the  first  time  period,  AT,  the  pressure  sensed  at  a  distance 
lQ  from  the  bubble  center,  B0,  is  higher  than  the  imposed 
one.  As  a  result,  a  second  fictitious  bubble  placed  at  this 
distance  from  B0  would  have  a  slower  growth  during  this  first 
period,  AT.  This  phenomena  is  however  reversed  in  the  second 
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phase  as  an  expansion  wave  is  generated  by  the  growing  bubble 
B0.  In  the  third  and  last  phase  a  compression  wave  increases 
the  driving  pressure  for  collapse  making  this  one  more  intense. 

In  Figure  4.10,  we  see  an  example  of  the  pressures 
generated  during  the  bubble  history  at  two  locations:  a)  the 
center  of  the  cloud  and  b)  the  center  of  one  bubble,  B* ,  in 
its  absense.  These  pressures  are  compared  with  those  generated 
during  the  growth  and  collapse  of  an  isolated  bubble  at  a 
distance  equal  to  the  spherical  cloud  radius.  We  have  selected 
the  case  of  a  finite  time,  AT  =  0.6,  pressure  drop.  The  same 
observations  made  while  interpreting  Figure  4.9,  can  be  repeated 
here  but  with  much  more  accentuated  values  of  the  pressure 
rises.  After  the  imposed  ambient  pressure  increases,  the 
nondimensional  pressures  generated  by  the  twelve  bubble  cloud 
are  first  positive,  then  a  pressure  expansion  period  is  observed 
for  1.9  <  t  <  3.4,  followed  by  a  high  pressure  surge  at  the  end 
of  the  collapse.  The  corresponding  bubble  radius  variation  with 
time  is  that  represented  in  Figure  4.8  (12  bubbles  AT  =  0.6). 

Figure  4.11  is  a  collection  of  the  results  obtained  from 
several  cases  studied.  The  maximum  nondimensional  pressures 
generated  during  the  cloud  collapse  are  represented  versus  the 
number  of  bubbles  in  the  cloud.  The  cumulative  effect  is 
obvious  since  the  values  obtained  vary  over  a  range  of  several 
orders  of  magnitude.  The  numbers  presented  should  not  be 
considered  accurate  since  other  scales  for  times,  pressure,  and 
lengths  are  needed  at  the  end  of  the  collapse.  Instead,  they 
are  presented  here  to  give  an  indication  of  how  tremendous 
pressures  can  be  generated  with  an  increasing  number  of 
interacting  bubbles,  and  to  give  an  idea  of  the  trend  of  this 
increase.  In  this  figure,  the  maximum  pressures  are  presented 
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at  the  cloud  center,  C,  at  B1  when  it  is  removed  and  at  a 

distance  rw  from  the  center  of  one  of  the  bubbles  in  the 
uo 

cloud . 

These  results  show  the  important  role  played  by  the  gas 
content  of  the  bubbles  which  was  neglected  in  March's  contribu¬ 
tions.  Increasing  PgQ  from  0.1  to  0.2  has  dramatically 
reduced  the  generated  pressures.  This  is  mainly  due  to  the 
cushioning  effect  of  the  gas  which  significantly  reduces  the 
velocities  attained  at  the  end  of  the  implosion. 

Another  very  interesting  observation  from  Figure  4.11  is 
that  the  imposed  pressure  variation  B  (pressure  drop  of  finite 
duration  followed  by  a  recompression)  moves  the  maximum 
pressures  generated  at  the  end  of  the  collapse  toward  much  lower 
values  than  for  pressure  jump  case  A.  This  effect  is  not  due  to 
the  apparent  higher  gas  content  in  this  case.  Indeed,  the  value 
of  Pg  to  consider  for  comparison  purposes  should  be  for  all 
cases  that  at  the  start  of  the  collapse  —  when  the  bubble  has 
its  maximum  volume.  For  example,  for  the  case  of  twelve  bubbles 
and  a  pressure  drop  ( AT  =  0.8,  Pgo  =  0.  53)  the  value  of 
Rmax/Ro  is  1*^3  (Figure  4.7).  Then,  accounting  for  the  gas 
expansion,  the  gas  pressure  at  the  beginning  of  the  collapse  is 
Pg  =  PgQ  (1.63)“4'2  =  0.07.  The  effective  gas  content 
is  thus  smaller,  and  since  the  value  of  e  is  bigger  (0.2  instead 
of  0.1),  the  observed  pressure  drop  is  intrinsically  related  to 
the  imposed  pressure  function.  The  pressure  attenuation 
observed  is  explained  by  the  initial  influence  of  the  cumulative 
effect  on  the  bubble  behavior  in  the  cloud  which  is  not  the  same 
at  the  start  of  the  growth  or  at  the  collapse  (Figures  4.3  and 
4.9).  In  the  classical  pressure  jump  case  the  initial 
cumulative  effect  is  to  prevent  the  bubble  size 
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from  being  small  when  the  collapse  pressure  surge  starts. 
Conversely,  the  initial  cumulative  effect  in  the  second  type  of 
imposed  ambient  pressure  (case  B)  is  to  reduce  the  size  of  the 
bubble  when  the  collapse  pressure  surge  occurs.  As  a  result  the 
potential  energy  for  collapse  in  case  B  is  relatively  smaller. 

Finally,  Figure  4.12  shows  the  influence  of  the  duration  of 
the  pressure  drop  on  the  maximum  pressures  generated.  By 
applying  the  same  reasoning  as  above  to  the  initial  gas 
pressures  one  can  conclude  that  the  increase  of  the  maximum 
pressure  with  AT  is  mainly  due  to  a  decrease  in  the  effective 
initial  gas  content  at  the  start  of  the  collapse  since  the 
maximum  bubble  radius  increases  with  AT. 

The  above  results  show  that  even  for  very  low  void 
fractions,  collective  bubble  collapse  can  generate  pressures 
orders  of  magnitude  higher  than  those  produced  by  single  bubble 
collapse.  This  would  tend  to  explain  the  observed  high  erosion 
intensities  and  the  bending  of  trailing  edges.  The  cumulative 
effect  comes  from  the  fact  that  the  interaction  increases  the 
driving  pressure  of  collapse  of  each  individual  bubble.  This 
augments  the  violence  of  its  implosion  and  thus  the  interaction 
with  the  other  bubbles.  Thus,  each  bubble  ends  its  collapse  not 
under  the  effect  of  a  pressure  of  the  same  order  as  the  ambient, 
but  orders  of  magnitude  higher.  This  cumulative  effect  would 
not  exist  if  the  void  fraction  is  high  enough  for  the  cloud  to 
behave  as  a  single  bubble.  This  leads  us  to  believe  that  there 
exists  a  critical  value  for  the  void  fraction  for  maximum 


erosion . 
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4 . 4  Comparison  to  Spark  Generated  Simulated  Multibubble 

Configurations 

We  describe  in  Section  8.0  some  experimental  tests 
conducted  to  simulate  symmetrical  bubble  cloud  configurations. 

A  spark-generated  bubble  method  is  used  in  combination  with  ! 

high-speed  photography  and  sound  measurements.  Bubbles  were  j 

generated  in  the  vicinity  of  complex  wall  configurations  such  1 

that  the  generated  bubble  and  its  images  in  the  walls  comprise  a  ' 

symmetrical  multibubble  set.  The  ratio  e  of  bubble  character¬ 
istic  size  to  interbubble  distance  was  varied  by  changing  the 
distance  between  the  electrodes  and  the  walls.  This  parameter 

I 

was  one  of  the  main  measurable  quantities  in  the  experiments. 

The  other  parameters  were  bubble  size  history  and  pressure  ; 

histories  at  selected  points  in  the  experimental  tank.  However,  j 

one  key  parameter  which  is  not  readily  measurable  is  the  gas 

pressure  inside  the  generated  bubble.  Since  knowledge  of 

Pgo  is  essential  for  any  quantitative  comparison  with 

experiment,  we  have  selected  to  use  it  as  a  free  parameter  and 

see  if  it  is  possible  to  match  the  experimental  results  by 

selecting  a  particular  value  of  PgQ. 

Figure  4.13  shows  an  example  of  the  influence  of  PgQ 
on  the  collapse  of  4  bubbles  (or  a  bubble  in  a  two-dimensional 
corner).  As  mentioned  earlier,  the  collapse  is  significantly 
softened  when  the  initial  gas  pressure  (or  the  initial  amount  of 
gas)  in  the  bubble  is  increased  (behavior  similar  to  the 
isolated  bubble  case).  The  bubble  oscillates  smoothly  for 
Pgo  0.15,  while  a  significant  violent  implosion  is 
observed  for  PgQ  <  0.10.  This  behavior  is  further  enhanced 
for  larger  values  of  e.  This  can  be  seen  in  Figure  4.14  for  the 
same  case  of  a  four-bubble  cloud  (PgQ  =  0.1),  and  in  more 
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detail  in  Figures  4.15  and  4.16  for  a  twelve  bubble 
configuration.  These  three  figures  show  the  behavior  of  a 
bubble  in  the  cloud  when  either  the  bubble  sizes  are  increased 
with  the  same  spacing  or  the  spacing  is  decreased  for  the  same 
bubble  size.  First,  the  classical  lenghtening  effect  of  the 
bubble  period  is  observed  (Figure  4.15,  e  =  0.01  and  Figure 
4.16,  e  =  0.02).  For  larger  values  of  e,  however,  the  same 
effect  as  an  increase  in  the  number  of  bubbles  (Figure  4.6)  is 
observed  —  namely  an  initial  slowing  down  of  the  collapse 
process  followed  by  a  very  significant  enhancement.  As  a 
result,  the  collapse  time  is  shortened.  For  larger  values  of  e, 
e  >  0.07b  in  Figure  4.15  and  e  >  0.14  in  Figure  4.16,  the 
asymptotic  method  obviously  fails.  A  reversal  of  the  collapse 
process  is  observed  and  cannot  be  easily  explained.  There  might 
be,  however,  some  correspondance  between  these  trends  and  the 
experiments  which  indicate  the  presence  of  an  optimum  value  of  e 
for  maximum  collapse  intensity  (see  Section  8.0). 

Figures  4.17  to  4.19  attempt  to  quantitatively  compare  the 
theoretical  model  and  the  spark-generated  bubble  experiments. 

As  mentioned  earlier,  the  difficulty  with  the  experiment  is  the 
inability  to  control,  or  even  to  simply  measure,  the  initial  gas 
pressure  in  the  bubble.  This  raises  some  questions  about  the 
accuracy  of  the  results  obtained.  With  this  in  mind  we  ran 
several  sets  of  numerical  cases  in  which  PgQ  and  the  number 
of  bubbles  was  maintained  constant  and  e  varied.  In  each  case, 
the  value  of  the  normalized  time  at  which  the  bubble  radius 
achieved  its  first  minimum  was  noted,  and  the  values  compared 
with  the  experimental  observations.  Figures  4.17  to  4.19  show 
plots  of  these  values  of  the  bubble  "first  period  of  oscilla¬ 
tion"  as  function  of  e.  Even  if  the  comparison  does  not  seem  to 
be  very  encouraging  at  first,  a  closer  look  enables  one  to  make 
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several  interesting  observations.  Oscillations  in  the 
theoretical  curves  indicating  strong  nonlinearities  appear  as 
for  the  experimental  curves.  This  is  particular ily  visible  in 
Figure  4.17  for  the  case  of  a  four-bubble  cloud.  These 
oscillations  are  par ticular ily  important  when  the  initial 
normalized  gas  pressure  PgQ  is  high.  The  combination  of  the 
experimental  points  and  the  theoretical  curves  seem  to  indicate 
that  for  the  four-bubble  case  PgQ  was  in  the  range  0.1  to 
0.2.  Examination  of  Figures  4.18  and  4.19  indicate  lower  values 
for  the  two  other  series  of  experiments  —  simulated 
eight-bubble  case  (bubble  in  a  three-dimensional  corner)  or 
twelve-bubble  case  (bubble  in  a  pyramidal  corner).  A  value 
lower  than  0.1  seems  to  apply  for  Figure  4.19,  and  PgQ  = 

0.05  seems  to  apply  to  Figure  4.18.  These  differences  with  the 
same  spark  generation  system  in  the  same  work  conditions  can  be 
related  to  different  degrees  of  water  degasing,  or  to  the 
rearrangement  of  the  electrodes  between  one  set  of  experiments 
and  another.  Scatter  and/or  oscillations  of  the  results  seem  to 
be  reduced  with  PgQ-  However,  the  general  trend  of  the 
existence  of  minimum  values  of  bubble  oscillation  period  or 
higher  collapse  pressure  (see  Section  8.0)  seems  to  be  confirmed 
by  both  theory  and  experiment.  Shima,  et  al .  (1983)  have  made 
similar  observations  for  another  range  of  values  of  e,  e  ^  1  . 

4 . 5  Study  of  Multibubble  Dynamics  in  a  Venturi  Cavitation 
Susceptibility  Meter 

4.5.1  Introduction 

The  development  of  techniques  to  measure  cavitation  nuclei 
distribution  in  laboratories  and  oceans  has  been  intensified  in 
recent  years.  Review  studies  on  the  subject  have  been  presented 
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recently  by  Shen  and  Peterson  (1983)  and  Billet  (1984).  One  of 
the  devices  being  developed  and  to  which  the  multibubble  study 
will  be  applied  here  is  the  Cavitation  Susceptibl ity  Meter 
( CSM) .  This  is  a  venturi  system  with  a  specially  tailored 
geometry  for  cavitation  nuclei  measurements  (Oldenziel,  1982  and 
Lecoffre  and  Bonnin,  1979).  Microbubbles  entering  the  venturi 
grow  at  the  throat  and  are  detected  either  optically  or 
acoustically.  In  the  first  case  they  have  to  exceed  a  certain 
size  to  be  perceived,  while  in  the  second  case  the  noise 
generated  during  the  implosion  of  these  bubbles  has  to  exceed  a 
certain  acoustical  level.  Therefore  in  this  method  a  "critical 
pressure",  Pcr ,  at  the  venturi  throat,  below  which  cavitation 
events  are  detected  by  the  CSM,  is  measured. 

The  classical  cavitation  number,  a,  (Thomas  number)  has 
been  widely  used  in  scaling  cavitation  inception  between  model 
and  prototype: 

P  -  P 

oo  V 

^  v  /  4  1  \ 


In  this  expression  P„  and  V„  are  the  characteristic  pressure 
and  velocity,  p  the  liquid  density,  and  Pv  the  liquid  vapor 
pressure.  This  definition  assumes  that  cavitation  inception 
occurs  when  the  liquid  pressure  drops  below  vapor  pressure. 
However,  it  has  been  recognized  for  a  long  time,  using  static 
equilibrium  theory,  that  there  exists  for  each  initial  nuclei 
size  a  critical  pressure,  usually  much  lower  than  pv,  below 
which  unstable  bubble  growth  occurs.  Since  in  a  liquid,  nuclei 
are  distributed  in  a  relatively  wide  range  of  sizes,  the 
definition  of  a  "liquid  critical  pressure"  can  only  be  done  by 
introducing  an  inception  criterion  such  as  a  threshold  amplitude 
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and  rate  of  the  detected  cavitation  "events."  It  has  been 
demonstrated  experimentally  that  if  py  is  replaced  by  such  a 
"critical  pressure"  of  the  liquid,  the  correlation  of  cavitation 
inception  observed  experimentally  between  different  model  sizes 
is  greatly  improved  (Keller,  1984).  It  is  implied  that  by  doing 
so  one  accounts  for  differences  in  the  nuclei  population  without 
having  to  determine  the  actual  bubble  size  distribution.  With 
the  CSM  a  "critical  pressure"  at  the  venturi  throat  is 
determined  by  varying  the  throat  velocity  and  counting  the 
occurrence  of  cavitation  bursting  optically  or  acoustically. 

If,  in  addition,  there  exists  a  thorough  theoretical  knowledge 
of  the  relationship  between  fluid  flow  characteristics ,  nuclei 
initial  size  and  needed  critical  pressures  to  initiate 
cavitation  events,  the  CSM  might  be  useful  as  a  bubble  nuclei 
size  detector.  Thi  section  attempts  to  improve  on  this 
knowledge  by  considering  bubble  dynamics  and  interactions 
between  bubbles  and  with  the  main  flow  into  the  venturi.  This 
investigation  is  to  be  added  to  earlier  studies  by  d'Agostino 
and  Acosta  (1983)  and  Shen  and  Gowing  (1984). 

Cavitation  susceptibility  meters  are  designed  in  such  a  way 
that  when  in  operation  there  is  only  one  bubble  in  the  venturi 
throat  at  a  given  time  (See,  for  instance,  the  paper  of 
d'Agostino  and  Acosta,  1983).  However,  in  practice  multibubbles 
appear  occasionally  in  the  venturi  and  invalidate  the  method  of 
detection.  In  this  section  multibubble  interaction  is  included 
in  the  study  of  cavitation  inception  in  the  CSM  venturi  system, 
and  its  influence  on  the  inception  criterion  is  investigated. 

Due  to  the  restriction  of  the  pipe  wall,  the  velocity  field 
and  the  pressure  field  in  a  venturi  system  are  altered  when  the 
microbubbles  begin  to  grow.  This  subject  has  been  studied  by 
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d'Agostino  and  Acosta  (1983)  for  a  single  spherical  bubble 
entering  from  a  reservoir  into  a  semi-infinite  pipe  of  constant 
cross-section.  In  the  present  paper  we  consider  the  case  where 
multiple  bubbles  are  present  in  a  venturi  and  we  derive  an 
expression  for  the  pressure  correction  due  to  the  wall 
restriction.  This  correction  is  included  as  an  option  in  the 
symmetrical  multibubble  code,  CLDMAIN1,  when  the  bubbles  are 
grouped  together  as  a  cloud.  The  growth  and  collapse  of  a 
single  bubble  and  multibubbles  in  a  venturi  system  with 
contraction,  throat,  and  diffuser  zones  of  various  diameters  are 
computed  considering  the  venturi  wall  restriction. 

4.5.2  Single  Bubble  Static  Equilibrium  and  Dynamics 

We  have  earlier  defined  the  critical  pressure  of  a  liquid 
as  the  pressure  (or  tension)  below  which  significant  bubble 
growth  or  activity  is  detectable.  In  fact,  as  can  be  deduced 
from  the  dynamic  stability  analysis  of  the  static  equilibrium  of 
isolated  spherical  bubbles,  there  is  a  different  critical 
pressure  for  each  initial  bubble  size  considered.  The  balance 
of  pressure  across  the  bubble  interface  can  be  written: 

P  =  Pv  +  Pg  -  2  y/R  ,  (4.2) 

where  P  is  the  ambient  pressure,  Pv  and  y  are  the  liquid  vapor 
pressure  and  surface  tension  and  R  is  the  equilibrium  bubble 
radius.  Pg  is  the  pressure  of  noncondensible  gas  inside  the 
bubble.  If  we  neglect  gas  diffusion  (see  Section  6)  we  can 
assume  that  the  gas  behaves  ideally  with  the  law  of  compression: 


Constant 


(4.3) 
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the  polytropic  constant  K  lies  between  1  (isothermal  case)  and 
cp/cv  (adiabatic  case).  The  relationship  between  the 
ambient  pressure,  P,  and  the  equilibrium  radius,  R,  is  therefore 
univocally  defined  for  a  given  reference  value  Rq,  or 
Pqo corresponding  to  a  reference  pressure  PQ: 


with 


»  =  P  - 

g  o 


P  + 
v 


2  Y 


(4.5) 


Figure  4.20  shows  an  example  of  the  static  equilibrium  curves 
P(R)  for  the  isothermal  case,  K  =  1. 


The  unstable  portions  of  these  curves  correspond  to 
dP/dR  >  0,  while  the  critical  pressure,  Pcr  and  radius,  Rcr , 
are  given  by  dP/dR  =  0.  These  critical  values  are  related  to 
each  other  and  to  the  initial  conditions  by  the  relations: 


R 


3K-1 

cr 
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Therefore  a  bubble  of  initial  size  Rq  at  the  ambient  pressure 
PQ  would  explosively  grow  (no  static  equilibrim)  if  the 
pressure  around  it  drops  below  Pcr.  One  should  note  that 
Pcr  is  always  lower  than  Pv  and  tends  toward  it  when  Rq  is 
very  large. 
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A  question  of  great  interest  to  the  venturi  CSM  as  an 
eventual  bubble  size  detector  is,  for  a  given  reference 
pressure  PQ  (depth  of  submergence)  and  throat  velocity  Vt, 
what  is  the  initial  radius  of  the  bubble  or  nucleus  which  would 
reach  the  critical  size  at  the  throat  pressure,  Pth?  All 
bubbles  of  larger  radius  would  then  be  unstable,  while  smaller 
bubbles  would  remain  undetectable.  In  terms  of  the  maximum 
throat  velocity  (accounting 
effects)  Pth  can  be  written 

Pth  Po  +  J  Cpmin 


for  pressure  loss  by  viscous 
P  Vfc2  ,  (4.8) 


and  the  relation  between  R 0  and  the  critical  throat  velocity, 
Vtc,  is 

pth  -  Pv  -  Po  +  2  0  Cpmin  Vtc  '  Pv 

-1 

2y(1-3K)  ^  C3KRo  (Po-Pv)]  3k_1 

3KR  3K  +  2y 

°  L  (4.9) 


For  a  given  critical  throat  velocity,  Vtc,  any  bubble  larger 
than  Ro  would  grow  explosively.  Conversely  for  a  given  bubble 
size  any  Vt  larger  than  the  critical  velocity  would  render  the 
pressure  at  the  throat  below  the  bubble  critical  pressure. 

Figure  4.21  shows  the  relationship  between  Rq  and  Vt  for 
various  ambient  pressures  and  for  K  =  1.4.  Note  the  sensitivity 
to  the  critical  velocity  for  relatively  "large"  radii  (e.g., 

Rq  >  1  um).  Very  small  changes  in  Vt  induce  large  changes 
in  Rq.  This  underlines  the  major  practical  problem  in  the  use 
of  the  CSM  as  a  bubble  size  detector.  Much  less  sensitivity 
would  be  obtained  if  Pth  were  the  variable  controlled  in  the 
experiments  (see  Figure  4.20). 
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The  original  concept  for  the  use  of  the  venturi  cavitation 
susceptibility  meter  was  based  on  the  above  static  equilibrium 
approach.  Actually,  the  above  reasoning  can  be  invalidated  when 
bubble  dynamics  is  taken  into  account  in  the  two  following 
configurations  (Darrozes  and  Chahine,  1983): 

a.  The  pressure  in  the  throat  drops  below  the  critical 
pressure  but  the  considered  bubble  does  not  remain  long  enough 
in  this  region  to  grow  explosively. 

b.  The  imposed  pressure  is  always  above  the  critical 
pressure  but  pressure  variations  are  great  enough  for  an  intense 
collapse  and  therefore  a  strong  acoustical  signal  to  be 
generated . 

In  order  to  account  for  the  above  phenomena,  the  complete 
bubble  dynamical  equation  must  be  solved  for  the  particular 
venturi  CSM  pressure  field.  In  the  case  of  a  single  bubble  of 
negligible  size  compared  to  the  throat  diameter,  the  classical 
Rayleigh-Plesset  equation  adequately  describes  the  bubble  radius 
variations.  The  external  pressure  in  this  equation,  P(t),  is 
the  pressure  "felt"  by  the  bubble  during  its  traverse  of  the 
venturi.  It  is  the  pressure  in  the  liquid  in  the  absence  of  the 
bubble  at  the  location  of  its  center.  The  bubble  position  x(t) 
is  obtained  from  the  bubble  translation  velocity  V^ft) 

x(t)  =  /£  Vb(t)  dt  .  (4.10) 

If  we  neglect  relative  velocity  between  the  bubble  and  the 
liquid,  Vb  is  the  fluid  velocity  which  varies  along  the 
venturi.  Actually  the  slip  velocity  could  be  important  and 
should  be  included  in  future  developments  of  this  study  (see 
Section  4.5.5). 
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When  the  bubble  growth  rate  is  large  enough  to  signifi¬ 
cantly  modify  the  liquid  flow  rate  into  the  venturi,  the 
pressure  field  in  the  venturi  has  to  be  modified  to  account  for 
the  resulting  pressure  variations.  Similarly,  the  Rayleigh 
Plesset  equation  alone  is  not  adequate  to  describe  the  case 
where  interaction  between  several  bubbles  crossing  simultane¬ 
ously  the  venturi  is  considered. 

4.5.3  Multibubble  Interaction  with  Venturi  Mean  Flow 

Under  the  basic  assumption  of  small  bubble  radius  size 
relative  to  bubble  spacing  and  to  flow  characteristic  length 
(for  instance  throat  radius),  in  first  approximation  each  bubble 
behaves  independently  and  the  basic  flow  is  undisturbed  by  the 
presence  of  the  cloud.  With  this  approach  and  at  this  order  of 
approximation,  the  pressure  field  to  which  each  bubble  in  the 
cloud  reacts  is  that  existing  at  the  location  of  its  center  in 
its  absence.  At  higher  orders  bubbles  interact  with  each  other, 
and  due  to  the  restricted  nature  of  the  flow  in  the  venturi  they 
modify  the  main  flow  into  the  venturi.  The  influence  of  the 
cloud  dynamics  on  the  basic  flow  becomes  important  when  the 
cloud  volume  rate  attains  a  significant  fraction  of  the  mean 
volume  flow  rate.  We  will  try  here  to  account  for  this 
correction  since  its  implementation  is  rather  simple. 

Let  r<  t)  be  the  total  volume  of  bubbles  present  at  time  t 
in  the  venturi  (these  bubbles  could  be  either  a  certain  number 
of  isolated  separated  bubbles  or  a  bubble  cloud).  The 
difference  between  the  upstream  and  downstream  flow  rates  into 
the  venturi  is  equal  to  r'-  With  a  unidimensional  approach  to 
the  venturi  flow  it  is  reasonable  to  write  that  the  mean  flow 
entering  the  venturi  is  then  modulated  as  follows: 
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V  =  a  -  a 
o  o 


r/Ac 


(4.11) 


where  Aq  is  the  area  of  the  venturi  at  its  entrance  and  UQ 
the  undisturbed  velocity  at  this  location.  a  is  the  fraction 
of  the  displaced  liquid  volume  rate  due  to  the  bubble  which 
propagates  upstream  during  the  bubble  growth.  a  needs  to  be 
determined  with  an  additional  equation  not  available  presently. 
We  will  assume  from  here  on  that  a  =  0.5.  If  all  bubbles  were 
concentrated  at  x  =  X  then,  between  the  venturi  entrance  (x  =  0) 
and  X,  the  velocity  in  the  venturi  can  be  related  to  the 
undisturbed  velocity  U(x)  by  an  relation  equivalent  to  (4.11), 
valid  in  the  small  perturbation  approach: 


V(x)  =  U ( x )  - 


y 

?A('x) 


X  <  X 


(4.12) 


If  in  addition  the  flow  is  assumed  to  be  potential,  then  the 
velocity  potential  can  be  written: 


<f>(x)  =  4o(x)  - 


x  r 

■“  2A(s) 


(4.13) 


Bernoulli's  theorem  between  x  and  a  field  point  upstream  where 
the  pressure  is  Pa  and  the  fluid  is  at  rest  can  be  written: 


_a  =  P(x) 
P  P 


34>(x) 

at 


1  Uuil 

2  I  ax 


(4.14) 


We  can  now  define  a  perturbation  pressure  p' (x)  such  that 


P(x)  =  ’st  (x)  +  p'  (x)  , 


(4.15) 


where  Pst  is  the  steady  state  pressure  in  absence  of  bubbles. 
Combining  (4.13),  (4.14),  and  (4.15)  one  obtains  the  following 
expression  for  p'(x)  where  a  higher  order  term  in  C*f// A)  2  has 
been  neglected. 


-  -•‘.V.--  .  -  L--  .  - 
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p' (x)  = 
p 


1L 


2A(  a) 


da  +  U(x) 


Jl 

2A(x) 


(4.16) 


In  a  matched  asymptotic  approach  the  imposed  pressure  field 
around  any  bubble  or  cloud  at  location  X  is  given  by  (4.15),  and 
for  the  corresponding  first  order  bubble  the  imposed  pressure 
term  combines  the  undisturbed  pressure  in  the  venturi  at  the 
bubble  location  and  the  correction  p'(x)  given  by  (4.16). 

If  the  bubble  cloud  is  composed  of  N  bubbles  of  equal  size 

then 


y'  =  4  ir  NR3/3 


(4.17) 


and  the  Rayleigh  Plesset  equation  becomes: 


p  [RR  +  ~  R2] 


+ 


il 

R 


+ 


P(X(  t)  )  +  2 tt  pN 


/R2RU(X)  ,  -X 
^  A(X)  + 


2RR  2+R  ^R 
A(  x ) 


(4.18) 


D'Agostino  and  Acosta  (1983)  derived  a  similar  equation  for  a 
semi-infinite  pipe  of  constant  cross  section  Aq.  Equation 
(4.18)  combined  with  the  system  of  Equations  (3.11)  can  be 
solved  numerically  using  a  multi-Runge-Kutta  procedure,  to 
investigate  interaction  of  bubbles  with  the  venturi  main  flow 
and  with  each  other. 


4.5.4  Numerical  Results  and  Interpretation 

A  series  of  computer  runs  were  made  using  the  above 
described  method.  The  objective  of  this  investigation  was  to 
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compare  the  behavior  of  an  isolated  spherical  bubble  with  that 
of  a  bubble  in  a  cloud.  Both  sets  of  cases  were  considered  for 
the  NSRDC  venturi  shown  in  Figure  4.22  with  and  without  pressure 
correction  due  to  the  modification  of  the  mean  flow  by  the 
presence  of  the  bubbles.  The  results  obtained  are  illustrated 
in  the  following  examples  and  figures.  Due  to  viscous  effects  a 
Cpmin  of  -1.22  was  measured  at  NSRDC  for  this  venturi  (see 
Shen,  et  al . ,  1984)  and  was  accounted  for  in  the  following 
computations . 

The  bubble  wall  motion  during  its  traverse  through  the 
venturi  is  illustrated  in  Figure  4.23  for  a  bubble  of  10  urn 
initial  radius.  In  most  cases  computations  were  started  at  the 
throat  entrance  (X  =  0.825  cm)  in  order  to  reduce  computation 
time.  This  is  acceptable  since  comparative  tests  have  shown 
that  no  significant  error  is  then  introduced.  A  typical  bubble 
radius  history  can  be  described  as  follows.  Near  the  throat 
entrance  the  bubble  radius  oscillates  for  a  certain  number  of 
cycles  which  decreases  with  increasing  flow  velocity.  Then,  if 
this  velocity  is  large  enough  a  continual  growth  is  observed 
until  the  bubble  reaches  the  venturi  expansion  area.  This  is 
followed  by  a  more  or  less  violent  bubble  implosion  which  may  be 
detected  acoustically.  Since  this  behavior  is  primarily 
controlled  by  a  balance  between  gas  pressure  inside  the 
bubble  and  the  imposed  ambient  pressure,  a  key  factor  is  the 
polytropic  constant,  K.  The  isothermal  case  (K  =  1 )  is  the 
most  easily  manageable  mathematically  and  has  therefore  been 
studied  by  many  investigators.  Its  use  is  justified  when  the 
bubble  wall  motion  is  slow  enough  so  that  the  gas  temperature 
remains  constant.  However;  the  value  of  K  can  vary  widely 
during  the  bubble  history  (see  Plesset  and  Prosperetti,  1977  and 
Section  6).  Figure  4.23  shows  the  importance  of  K  by  comparing 
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the  results  for  the  two  classical  values  of  K:  1  and  1.4. 

Larger  maximum  radii  are  obtained  with  K  =  1.  Two  velocities  at 
the  throat  are  also  considered  in  Figure  4.23.  Here  we  see  that 
for  a  slight  variation  of  the  maximum  bubble  radius  can  be 
doubled.  Therefore,  if  the  criterion  of  cavitation  inception  is 
thatbubble  radius  exceeds,  let  us  say,  150  urn  then  Vt  =  25.0 
m/s  would  be  below  inception  and  Vt  =  25.2  m/s  would  be 
slightly  above  it.  The  static  equilibrium  curves  (Figure  4.21) 
derived  from  Equation  (4.9)  give,  however,  a  lower  value  for  the 
critical  throat  velocity,  V^.  =  24.95  m/s.  This  relatively 
small  difference  of  about  0.8  percent  is  significant  for  the 
precision  needed  in  the  measurement  of  Vt  in  order  to  use  the 
CSM  venturi  effectively.  Note  that  if  the  throat  velocity  at 
inception  was  24.95  m/s  the  critical  pressure  would  be  -500 
pascals  while  one  obtains  -8000  pascals  with  25.2  m/s.  The 
difference  in  Pcr  is  therefore  very  significant  between  the 
static  and  dynamic  approaches. 

Figures  4.24  and  4.25  show  the  maximum  size  attained  by  a 
bubble  crossing  the  venturi  CSM  as  a  function  of  the  velocity  at 
the  throat.  Bubbles  of  different  initial  sizes  at  the  same 
initial  reference  pressure  are  considered.  Single  bubble 
results  are  compared  with  symmetrical  multibubble  configura¬ 
tions.  Several  very  interesting  observations  can  be  made  from 
studying  these  figures  and  the  corresponding  numerical  results. 
Bubble  size  scaling  effects  are  observed  as  follows.  At  the 
lower  velocities  bubbles  exhibit  an  oscillatory  behavior  all  the 
way  through  the  venturi  (e.g.,  Vt  <  25.1  m/s  for  Rq  =  5  urn; 

Vt  <  24.8  m/s  for  Rq=  20ym).  Their  radius  oscillates 
between  a  fraction  of  Rq  and  a  maximum  radius  of  a  few  times 
Rq.  At  the  higher  velocities  the  bubbles  experience  an 
explosive  growth  followed  by  a  strong  implosion.  Bubble  size 
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scaling  effects  are  expressed  by  the  fact  that  while  initially 
larger  bubbles  have  a  larger  growth  rate  at  the  low  venturi 
throat  velocities,  the  trend  is  reversed  at  higher  velocities 
(Figure  4.24).  The  result  of  this  scaling  effect  is  that  at 
large  velocities  (e.g.,  =  25.8  m/s  for  the  cases 

considered)  single  isolated  bubbles  of  5,  20,  and  50  ym  initial 
size  attain  the  same  maximum  size.  This  could  be  seen  on  Figure 
4.24  where  the  dimensional  variables  are  represented  for  the 
same  data  points  as  in  Figure  4.25.  This  result  is  also  related 
to  the  fact  that  the  considered  venturi  has  a  long  enough  throat 
to  allow  the  above  bubbles  to  grow  to  significant  sizes. 

The  effects  of  bubble  interactions  and  restricted  mean  flow 
correction  show  up  in  a  similar  way  when  single  bubble  and 
multibubble  results  are  compared.  At  lower  velocities  bubble 
oscillations  in  the  throat  area  close  to  the  entrance  are 
amplified  when  the  collective  behavior  of  the  bubble  cloud  is 
considered  (see  Figure  4.27  for  details).  As  a  result  larger 
maximum  sizes  are  obtained  with  multibubbles  but  also  very  early 
collapses  are  achieved  in  the  throat  section  before  reaching  the 
expansion.  At  higher  velocities  the  opposite  phenomenon  is 
observed.  An  explosive  bubble  growth  is  then  achieved,  and  the 
multibubble  effect  is  to  inhibit  the  growth  rate  of  each  bubble 
in  the  cloud.  This  is  mainly  due  to  a  reduction  of  the  pressure 
drop  around  each  bubble  induced  by  the  motion  of  the  other 
bubbles.  One  can  notice  in  both  Figures  4.24  and  4.25  that 
while  multibubble  effects  weakly  modify  the  5  micron  bubble 
curves,  a  dramatic  change  is  seen  for  the  20  micron  bubble.  In 
the  latter  case  bubble  interactions  are  much  stronger  and  shift 
the  behavior  of  the  bubble  from  an  explosive  growth  to  strong 
oscillations  which  induce  an  early  collapse  even  for  velocities 
as  high  as  25.6  m/s.  At  higher  values  of  t  (bubbles  closer  to 
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each  other)  this  effect  can  be  seen  even  for  the  5  ym  bubble  at 
Vt  =  25.8  m/s  (case  e  =  0.05  on  Figures  4.24  and  4.25). 

Figures  4.24  and  4.25  could  also  be  used  to  determine 
dynamic  cavitation  inception.  This  could  be  done  by  using 
either  an  "acoustical"  criterion,  %ax/Ro  greater  than  a 
certain  value,  or  an  "optical"  criterion,  the  actual  bubble  size 
exceeding  a  given  value.  If  we  adopt  criteria  used  earlier  in 
the  literature  we  have  the  following  value  for  at  inception 
as  compared  to  static  equilibrium  values. 


Criterium 

Ro 

p  m 

Static 

m/s 

Vfc  Single  Bubble 
m/s 

V  5-Bubble 

m/s 

R/R  >  1 0 
o 

5 

25.08 

25.17 

25.13 

R  >  1  50  um 

5 

25.08 

25.24 

25.22 

R/R  >  10 
o 

20 

24.89 

25.13 

25.70 

R  >  150um 

20 

24.89 

25.08 

25.62 

From  the  above  table  one  can  notice  that  predictions  with 
static  theory  are  conservative.  For  a  given  bubble  size  they 
show  earlier  inception.  Similarly,  for  a  given  velocity,  static 
theory  predicts  the  activation  of  a  broader  range  of  nuclei. 

For  instance,  at  25.08  m/s  all  bubbles  larger  than  5  ym  are 
predicted  to  be  unstable  while  dynamics  show  that  only  those 
above  20 ym  will  be  observed  with  the  optical  criterion.  Larger 
differences  between  statics  and  dynamics  result  for  the  case  of 
a  multibubble  system.  One  should  however  recognize  that 
differences  in  Vt  are  relatively  small  especially  when 
considering  present  practical  considerations  in  the  control  of 
this  velocity.  On  the  other  hand,  note  that  the  resulting 
critical  pressure  varies  in  a  very  wide  range  between  the 
various  approaches.  The  difference  between  these  values  is  very 
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import  ant  to  account  for  in  cavitation  tunnel  tests  where  the 
ambient  pressure  is  low.  This  encourages  effort  toward  a  better 
control  of  or  Pcr.  This,  as  well  as  additional 
information  given  below,  might  explain  the  discrepencies  between 
optical  and  CSM  measurements  of  bubble  distributions.  Large 
size  bubbles  in  a  cloud  formation  are  not  excited  at  the 
velocities  expected  from  static  theory  and  are  therefore  not 
counted  when  interpreting  experimental  CSM  results. 

Figures  4.26  and  4.27  illustrate  some  details  of  the  results 
presented  in  Figures  4.24  and  4.25.  The  influence  of  the 
velocity  at  the  venturi  throat  on  the  dynamics  of  a  bubble  of 
initial  radius  5  urn  can  be  seen  in  Figure  4.26.  One  can  observe 
an  earlier  explosive  growth  of  the  bubble  at  the  higher  velocity 
as  well  as  the  achievement  of  a  greater  maximum  bubble  radius. 
The  bubble  reaches  its  maximum  size  in  the  venturi  diffuser  zone 
after  leaving  the  throat.  This  is  due  to  its  response  time  to  a 
change  in  the  ambient  pressure.  At  higher  velocities  the 
maximum  is  achieved  farther  downstream.  The  influence  of 
multibubble  interaction  and  of  restricted  mean  flow  modification 
on  the  bubble  dynamics  can  also  be  seen  in  Figure  4.26.  The 
explosive  growth,  as  observed  in  Figures  4.24  and  4.25,  is 
inhibited  by  multibubble  interaction  more  significantly  at  the 
higher  velocities.  The  location,  x,  of  maximum  radius  and  the 
location  of  the  collapse  are  then  closer  to  the  throat  exit. 

This  effect  increases  with  the  ratio  between  bubble  size 
and  distance,  e.  For  a  larger  t  the  influence  of  oscillations 
is  very  much  increased  and  the  bubble  collapses  much  earlier  in 
the  venturi  throat  (e.g.,  for  e  =  0.05,  the  collapse  is  x  =  1.30 
cm  instead  of  x  =  2.8  cm). 
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Figure  4.27  shows  the  influence  of  the  initial  size  of  a 
bubble  on  its  dynamics  and  the  modification  of  the  behavior  by 
interbubble  and  mean  flow  interactions.  All  cases  are 
considered  for  the  same  ambient  pressure  and  throat  velocity. 

The  larger  bubbles  achieve  their  maximum  size  further  downstream 
inside  the  venturi  diffusion  section.  When  interactions  are 
considered  the  bubble  of  initial  radius  Rq  =  5  um  has  its 
behavior  moderately  modified,  similar  to  what  is  shown  in  Figure 
4.26.  However,  for  the  20  um  bubble  the  modification  of  the 
behavior  is  dramatic.  Large  oscillations  are  induced  through 
interactions  and  an  early  collapse  at  X  =  1.35  cm  occurs  inside 
the  venturi  throat.  Such  a  bubble  cloud  would  not  be  detected 
optically  and  probably  not  acoustically  if  the  intensity  of  the 
collapse  is  not  strong  enough. 

For  the  present  user  of  the  Venturi  Cavitation 
Susceptibility  Meter,  the  curves  relating  the  initial  radius  of 
a  detectable  bubble  to  the  critical  velocity  at  the  throat  are 
the  most  useful.  Figure  4.21  showed  those  curves  based  on 
static  equilibrium  considerations.  Any  bubble  of  initial  radius 
larger  than  that  given  by  these  static  predictions  would  grow 
explosively  at  the  corresponding  velocity.  The  discussions 
above  pertaining  to  Figures  4.24  to  4.27  have  shown  that  the 
problem  is  more  complex  when  dynamics  and  interactions  are  taken 
into  account.  This  fact  can  be  illustrated,  as  in  Figure  4.28, 
by  comparing  Rq  versus  Vt  obtained  by  both  static  and 
dynamic  considerations.  Here  one  can  observe  significant 
discrepancies  at  the  larger  bubble  radii  and  also  at  the  higher 
velocities.  Again  these  theoretical  differences  are  relatively 
small  when  expressed  in  terms  of  Vt  but  reflect  very  large 
differences  in  the  critical  pressures.  The  most  interesting 
results  seems  to  be  the  presence  of  a  minimum  in  the  dynamic 
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curve.  This  implies  that  at  a  given  velocity  above  a  minimal 
value  (25.1  m/s  in  the  case  of  Figure  4.28)  only  a  finite  range 
of  bubble  radii  are  excited.  Compare  this  prediction  with  the 
conclusion  from  statics  that  all  bubbles  above  a  critical  radius 
become  active.  This  may  also  explain  the  experimentally 
observed  lower  number  of  detected  bubbles  in  the  venturi  CSM  as 
opposed  to  the  scattering  method. 

The  information  obtained  in  this  work  on  bubble  size 
history  in  the  venturi  as  a  function  of  the  imposed  flow 
conditions  and  initial  bubble  radii  should  be  complemented  and 
could  be  used  to  determine  nuclei  population.  This  could  be 
implemented  if  the  venturi  CSM  is  instrumented  to  measure  bubble 
sizes  in  the  venturi  at  several  locations.  This  could  be  done 
by  simple  optical  sensors. 

The  analyses  should  be  refined  and  extended  to  include 
bubble-liquid  slip  velocity  and  a  fine  description  of  bubble 
collapse.  Slip  velocity  could  modify  the  bubble  time  response 
and  quantitatively  change  the  results  obtained.  A  precise 
description  of  the  implosion  would  allow  us  to  accurately 
compute  the  pressures  generated  and  therefore  would  be  very 
helpful  for  acoustic  detection. 

4.5.5  Extension  of  Study  to  the  Flow  Around  a  Hydrofoil 

The  program  developed  to  study  bubble  behavior  in  a 
"slowly-varying"  pressure  field  (relative  to  the  bubble  period 
of  oscillation)  has  been  extended  from  the  venturi  case  to  that 
of  a  hydrofoil.  The  program  is  CLDMAIN2.  The  user  inputs  the 
pressure  profile  in  terms  of  coefficients  of  a  polynomial.  The 
study  differs  from  that  of  the  venturi  by  the  fact  that  gravity 
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forces  have  been  included.  Therefore,  slip  velocity  between  the 
bubble  and  the  surrounding  liquid  had  to  be  included.  This  is 
obviously  necessary  if  the  bubble  dynamics  in  the  wake  are  to  be 
monitored.  The  presence  of  gravity  introduces  also  a 
modification  of  the  input  pressure  field.  The  pressure  driving 
the  bubble  behavior  is  now  dependent  on  the  position  of  the 
bubble  relative  not  only  to  the  foil  but  also  to  a  free  surface. 

The  balance  of  forces  on  the  bubble  is  at  each  instant 
between  buoyancy,  drag,  and  acceleration  forces.  This  can  be 
written  for  a  bubble  of  center  B(t),  radius  a(t)  and  velocity 
VB(t): 

dV 

§  TT  a3  (  P  jL  +  2Pg)  =  4jrMa(VB-VJl)-2ira3Vp  + 

2  ir  a2  a*  (Vt  -  VB)  .  (4-19) 

Here  u  is  the  liquid  viscosity,  and  pg  are  the  liquid  and 
gas  densities,  Vp  is  the  pressure  gradient  and  V£  is  the 
liquid  velocity.  These  two  are  related  by  the  approximate 
relation 

Vp  =  dV  £/dt  ,  (4-20) 

and  for  the  present  case  (bubble  behind  foil) 

Vp  =  -  pg  .  (4-21  ) 


Therefore,  Equation  (4.19)  can  be  written  in  terms  of  the 
relative  velocity,  Vt. 
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Vt  =  -  [6  -2  +  3  f  ]  Vt  +  2  g  ,  (4-22) 

cl 

where  v  is  the  liquid  kinematic  viscosity  and 

vt  =  VB  "  vi  ♦  (4-23) 

Equation  (4.22)  was  solved  using  a  Runge-Kutta  procedure  similar 
to  those  used  for  the  various  bubble  radius  components  and  was 
implemented  in  CLDMAIN2.  This  implementation  is  in  fact  very 
similar  to  that  needed  for  the  general  case  of  a  slip  velocity 
(due  to  inertial  and  drag  effects)  between  the  bubble  and  the 
surrounding  liquid. 

4 • 6  Asymmetric  Bubble  Cloud  Configurations 

In  this  section  we  present  a  few  illustrations  of 
multibubble  dynamics  when  all  bubbles  in  the  cloud  do  not  behave 
identically  due  to  symmetry.  To  do  so  we  have  used  the 
developed  code  MULTIBBL  (see  flow  chart.  Figure  4.2)  in  which 
the  position  and  size  of  the  bubbles  are  entered  using  cartesian 
or  spherical  coordinates.  Figures  4.29  and  4.30  show  the 
results  for  the  simple  example  of  a  configuration  of  three 
bubbles  of  different  sizes.  The  comparison  between  the  bubble 
behavior  when  isolated  and  in  the  presence  of  the  two  other 
bubbles  shows  interesting  results.  When  two  large  bubbles 
collapse  in  the  presence  of  a  smaller  one,  the  rate  at  which 
their  radii  decrease  with  time  is  slightly  reduced  (Figure  4.29) 
and  the  oscillation  period  is  increased.  However,  the  behavior 
of  the  smaller  bubble  is  dramatically  modified.  Following  the 
end  of  the  first  collapse,  the  rebound  of  the  bubble  (due  to  the 
presence  of  noncondensibles)  is  much  stronger  and  is  presumably 
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intensified  by  the  perturbation  pressure  field  due  to  the  two  ! 

other  bubbles'  implosion  (see  Figure  4.4).  Later  the  small  J 

i 

bubble  has  a  weak  growth  and  implodes  strongly  for  t-0.93  when 
the  local  pressure  field  attains  a  peak. 

The  case  of  two  small  bubbles  and  one  large  one  (Figure 
4.30)  shows  the  same  trends.  However,  in  this  case  the  large  j 

bubble  is  much  less  influenced  by  the  presence  of  two 
smaller  ones  than  by  the  presence  of  one  smaller  and  one 
identical  bubble  as  in  Figure  4.29.  The  deviation  of  the  small 
bubbles'  radii  from  the  isolated  case  is  also  much  less 
dramatic . 

Figures  4.31  and  4.32  address  the  question  of  the  influence 
of  cloud  geometry  (bubble  space  distribution)  on  cloud 
dynamics.  Three  multibubble  configurations,  all  having  the  same 
global  void  fraction  (total  bubble  volume  over  volume  in  which 
they  are  distributed),  are  compared.  In  the  three  cases,  twelve 
bubbles  of  equal  initial  sizes  (R0=1)  are  located  on  or  inside  a 
sphere  of  nondimensional  radius  15.  In  the  first  case,  the 
centers  of  all  twelve  bubbles  are  located  on  the  sphere  of 
radius  15.  In  the  second  case,  six  are  centered  on  the  same 
sphere  while  the  six  others  are  on  a  concentric  sphere  of  radius 
7.5.  The  third  configuration  is  composed  of  four  bubbles 
centered  on  the  outer  sphere  of  radius  15,  four  others  are 
centered  on  a  sphere  of  radius  10,  and  the  innermost  layer  of 
four  bubbles  is  on  a  sphere  of  radius  5.  Figures  4.31  and  4.32 
show  that  the  cumulative  effect  is  greatly  enhanced  for  the 
inner  layer  of  bubbles.  The  initial  weakening  of  the  implosion 
as  well  as  its  later  amplification  are  more  pronounced  (even  for 
the  outer  shell  of  bubbles)  when,  for  the  same  global  void 
fraction,  the  number  of  shells  is  increased.  Similarly  in  the 
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same  cloud  the  effects  are  more  pronounced  for  the  inner-most 
layer  of  bubbles.  Apparently,  the  predominant  factor  is  the 
relative  bubble  size  and  spacing  (parameter  e=rb0/Jl0  in  the 
theory) .  In  both  figures  the  formation  of  a  reentering  jet  at 
each  bubble’s  wall  occurs  when  the  radius  becomes  zero.  The 
negative  values  seen  in  the  figure  are  meaningless,  and  the 
computation  should  be  stopped  at  that  point.  One  can  notice 
that  this  occurs  much  earlier  for  the  shell  closest  to  the  cloud 
center  (Figure  4.32). 

Figure  4.33  is  intended  to  show  the  capabilities  of  the 
program  MULTIBBL  as  well  as  the  complexity  of  the  bubble 
behavior  for  a  random  bubble  configuration.  Here,  twelve 
bubbles  have  been  located  randomly  in  a  space  enclosed  within  a 
sphere  of  radius  7.  Similarly,  the  initial  bubble  radii  were 
chosen  in  a  random  fashion  between  0  and  1 .  Both  these  are 
shown  via  the  table  in  Figure  4.34.  We  will  not  try  to  dwell 
too  much  on  the  interpretation  of  this  figure  since  the 
randomness  makes  any  rapid  analysis  empirical.  As  a  general 
statement,  we  can  say  that  while  each  bubble  has  the  tendency  to 
collapse  at  its  own  period  of  oscillation,  the  interactions  make 
the  behaviors  at  the  rebounds,  and  the  following  collapses  quite 
independent  of  the  initial  radius.  In  addition,  the  smaller  the 
initial  bubble  size,  the  more  influence  of  the  interaction  can 
be  seen. 

Figures  4.35-4.37  consider  the  case  of  the  growth  and 
collapse  of  a  symmetrical  configuration  of  six  bubbles  near  a 
solid  wall.  A  pressure  drop  of  finite  duration,  AT  =  0.8,  has 
been  imposed  on  the  six-bubble  cloud.  Four  of  the  bubbles  are 
in  a  plane  parallel  to  the  wall  and  at  a  normalized  distance  of 
12.5  from  this  wall.  The  two  last  bubbles  are  on  a  line 
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perpendicular  to  the  wall  and  are  centered  respectively  at  5  and 
20  bubble  radii  from  the  wall.  The  presence  of  the  wall  creates 
the  asymmetry  and  shows  three  different  behaviors  of  the 
initially  identical  bubbles.  The  geometrical  configuration 
being  such  as  described  above,  the  four  bubbles  in  the  plane 
parallel  to  the  solid  wall  behave  identically  (bubble  3  in  the 
figure).  The  two  remaining  bubbles  on  the  axis  perpendicular  to 
the  solid  wall  have  different  behaviors.  The  bubble  closest  to 
the  wall  is  bubble  No.  t  in  the  figures  and  the  farthest  away 
is  bubble  No.  6  in  the  figures.  Figures  4.35  and  4.36  indicate 
that  the  four  bubbles  in  the  midplane  seem  to  be  the  most 
influenced  by  bubble  interactions.  They  are  followed  by  the 
bubble  closest  to  the  wall.  This  result  is  probably  due  to  the 
fact  that  the  distance  to  the  wall  is  not  small  enough  from  the 
image  cloud  to  influence  the  closest  bubble  in  the  analyzed  case 
(bubble  radii=1 ,  location  of  centers  on  sphere  of  radius=7.5, 
and  distance  from  wall=5). 

The  comparison  between  Figures  4.35  and  4.36  shows  the 
influence  of  a  finite  sound  speed  in  the  liquid  on  the  bubble 
radius  history.  For  a  normalized  sound  speed  of  30  (actual 
speed  =  30  1 0/  (  r^  Q  yj  p/  AP )')  ,  Figure  4.36  shows  a  weakened 
interaction  compared  to  the  infinite  sound  speed  case.  As 
discussed  in  Section  4.2,  this  is  due  to  the  fact  that  a  finite 
sound  speed  induces  a  time  delay  between  the  emission  of  the 
pressure  oscillations  from  one  bubble  and  its  arrival  to 
another.  The  flow  due  to  a  point  source  of  intensity  q(t),  for 
instance,  at  a  distance  r  is  r"2q(t-r/c)  instead  of 
r-2q(t).  The  result  shown  in  Figures  4.35  and  4.36  is, 
however,  not  general  since  a  time  delay  could,  depending  on  its 
amount  and  the  shape  of  the  various  bubble  pressure  functions, 
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act  in  such  a  way  as  either  to  intensify  or  dampen  the 
implosion . 

Figure  4.37  shows  the  shape  of  the  bubbles  at  different 
times  during  the  implosion  process.  The  three  bubble  behavior 
types  described  above  (four  bubbles  in  the  middle  plane, 
farthest,  and  closest  bubble  to  the  wall)  are  shown  side  by 
side.  These  profiles  give  a  better  idea  of  the  behavior  and 
deformation  of  the  three  bubble  types. 
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5.0  INCLUSION  OF  HEAT  TRANSFER  AT  THE 
BUBBLE  WALL:  SUPERHEATED  LIQUIDS 

5.1  Introduction 


Many  modern  processes  deal  with  various  fluids  in 
conditions  where  both  heat  transfer  effects  and  inertia 
contribute  in  controlling  the  bubble  behavior.  Examples  of  such 
fluids  are  hydrocarbons,  liquid  metals,  cryogenic  fluids, and 
demineralized  hot  water  at  temperatures  as  high  as  300 °C.  Heat 
transfer  boiling  or  cavitation  appears  with  these  liquids  in 
such  applications  as  high  speed  flows  of  sodium-cooled 
fast-breeder  reactors  in  nuclear  power  engineering,  circulation 
of  cryogenic  liquid  in  pumps  in  aerospace  engineering,  and  flow 
of  hot  water  in  nozzles  and  tubes  in  steam  power  plants. 
Accidents,  such  as  loss  of  vacuum  insulation  in  cryogenic 
storage  tanks  and  loss  of  coolant  in  nuclear  power  plants,  are 
sources  of  boiling  nucleation  and  of  major  safety  concern 
(Plesset,  1980). 

In  the  previous  sections,  we  investigated  analytically  and 
numerically  the  collapse  of  a  bubble  cloud  due  to  an  increase  of 
the  ambient  pressure  neglecting  heat  transfer.  A  cumulative 
effect  was  shown  leading  to  pressures  generated  during  the 
collapse  significantly  larger  than  would  be  computed  by  adding 
the  effects  of  individual  bubbles.  This  explained  the 
observations  of  the  bent  trailing  edge  of  propellers  subjected 
to  cloud  cavitation.  In  this  section,  we  extend  the  singular 
perturbation  approach  earlier  developed  to  study  the  cases  where 
heat  transfer  effects  cannot  be  neglected.  We  then  investigate 
numerically  the  growth  of  a  bubble  cloud  in  a  superheated  fluid 
following  a  sudden  depressurization.  Both  a  general  approach 
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and  a  thermal  boundary  layer  approximation  are  studied 
analytically,  and  methods  of  numerical  solution  are  described. 
Numerical  computations  are  then  conducted  only  for  a  symmetrical 
cloud  bubble  configuration  with  the  boundary  layer  a^oroximation 
and  when  deviations  from  sphericity  are  very  moderate. 

The  general  equations  of  the  problem  have  been  presented  in 
Section  2.0  and  will  not  be  rewritten  here.  Similarly,  we 
consider  a  singular  perturbation  theory  as  in  Section  3.0 
assuming  a  small  ratio  of  bubble  size  to  bubble  interdistance. 

5 . 2  Singular  Perturbaton  Approach 

5.2.1  Normalizations 


In  order  to  make  asymptotic  expansions  (and  thus  to  compare 

orders  of  magnitudes)  an  accurate  choice  of  characteristic  scale 

variables  is  fundamental.  For  the  length  scales,  we  chose  as  in 

Section  3.0  the  bubble  characteristic  radius  rw  in  the 

uo 

inner  problem,  and  the  interbubble  distance  in  the  outer. 

However,  the  relationship  between  rfc>0  and  the  characteristic 
initial  bubble  radius,  Rq,  is  not  obvious.  Indeed,  while  in 

the  case  of  bubble  collapse,  the  bubble  radius  stays  of  order 
Rq  in  the  mathematical  sense,  (R( e)  =  0(Ro)  if  there  exists 
a  constant  A  independent  of  e  such  that  J  R |  <  X |  Rq  J ) ;  such  is 
not  the  case  for  the  bubble  cloud  growth  studied  here. 

Therefore  rbQ  is  chosen  arbitrarily  much  larger  than  Rq 
but  such  that  the  inequality 

r.  /I  =  e  <<  1  (5.1) 
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is  valid.  Consequently,  the  results  of  the  computations  will  be 
valid  only  as  long  as  the  radius  of  any  bubble  in  the  cloud  does 
not  greatly  exceed  r^. 

The  time  scale  choice  is  simple  once  rbQ  is  known.  In 
the  case  of  a  significant  pressure  drop,  as  for  the  problem  of 
sudden  depressurization  of  a  superheated  system,  this  time  scale 
is  related  to  the  pressure  drop  AP,  through 

tq  =  rb  y  p/ AP1  .  (5.2) 

°  o 


AP  could  also  be  the  order  of  magnitude  of  the  imposed  pressure 
fluctuations  when  P«,(t)  is  a  prescribed  function  of  time. 


As  mentioned  earlier,  in  both  the  "inner"  and  the  "outer" 
regions,  the  flow  in  the  first  approximation  is  that  due  to  a 
distribution  of  dynamic  sources  and  heat  sinks.  The  character¬ 
istic  strength  of  the  dynamic  sources  is  qQ  =  ^^q3/t0, 
and  depending  on  whether  one  considers  the  "inner"  or  the 
"outer"  problem,  the  resulting  velocity  potential,  $,  has  the 
scales : 


r 


iOUt 


=  r 


.  /  l  T 

b  o  o 
o 


(5.3) 


Since  the  maximum  temperature  drop  occurs  near  the  bubble 
wall, and  since  a  lower  bound  for  this  temperature  is  the  boiling 
temperature  of  the  liquid,  T^ ,  at  the  imposed  ambient 
pressure,  P„,  the  temperature  departure  from  T*,  is  scaled 
with  the  amount  of  superheat,  (T*  -  T^) . 


k 
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Wit  h  these  character istic  scales,  nondimensional  variables 
all  of  order  unity  are  introduced  through  the  following 
definitions,  where  bars  denote  outer  nondimensional  variables 
and  tildes  inner  ones: 

r  =  r/rb  ;  r  =  r/  tQ 
o 

,  ,  .  .  • 

l1J(t)  =  b(t)  =  bT0/r.  ; 

D  0 

(5.4) 

t  =  I  =  t/  T o ;  T  =  T  =  T/O^-Th); 

p(  t)  =  p(  t)  =  p(  t)/ AP  ; 

*(t)  =  -fin(t)/4on  ;  *(t)  =  *°Ut(t)/4°Ut. 

Each  of  the  unknowns,  X,  is  then  expanded  in  a  power  series  of  e 
as  follows: 

X  =  X  0+  eX  i  +  e  2X  2  +  e  3X  3  +  0(e3).  (5.6) 

5.2.2  First  Order  of  Approximations  ( e°) 

When  e  goes  to  zero,  the  distance  between  bubbles  goes  to 
infinity,  interactions  vanish,  and  in  the  absence  of  a 
slip  velocity  between  the  test  bubble  and  the  surrounding  fluid, 
the  only  boundary  condition  at  infinity  is  the  imposed  ambient 
pressure  variation  P„(t).  The  "inner  problem"  is  therefore 
spherically  symmetrical  and  its  solution  is  given  by  the  well- 
known  Rayle igh-Plesset  equation.  This  can  be  written  with  the 
superscript  (i)  omitted  for  convenience: 
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a  a  + 
o  o 


i  =  -pjt)  +  it  (t)  - 


2W 


(t)  -  . 


(5.7) 


a  o 


The  nond imens ional  parameters  are  defined  by  the  relations 


P  (t)  = 


Tro(t)  = 


P  (t)  -  P  Jo) 
_ 


pv(t)  -  Pv(o) 
AP 


PJ o)  -  Pv(o) 


AP 


#  W. 


-1  _  Y(o) 


(5.8) 


r  b  ^ 
o 


y(  t)  =  y( t)/ y(o) . 


Here  Y(t)  and  Pv(t)  are,  respectively,  the  surface  tension 
coefficient  and  the  liquid  vapor  pressure  at  the  bubble  wall 
temperature  at  time  t.  The  initial  equilibrium  condition  at  the 
bubble  interface  is 

+  2We-1/H0  =  0  (5.9) 


For  a  given  P^tt),  equation  (5.7)  can  be  solved  for  the 
variations  of  the  bubble  radius,  a01(t).  This  allows  the 
subsequent  determination  of  the  higher  order  approximations  of 
the  bubble  radius. 


When  the  temperature  at  the  surface  of  the  bubble  departs 
significantly  from  the  ambient  temperature,  it  is  necessary  to 
couple  equation  (5.7)  with  the  energy  equation  to  account  for 
the  dependence  of  py  and  y  on  temperature.  At  this  order,  the 
problem  is  spherically  symmetrical,  and  the  energy  equation, 
(2.5),  reduces  to  the  following  nond imensional  equation: 


a^a  3T 
o  o 


(5.10) 
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where  the  Peclet  number,  Pe ,  is  the  ratio  of  the  thermal  dif¬ 


fusion  time,  rb0Z/D,  to  the  bubble  characteristic  time  tqi 


"  rb  /(D  V 

o 


(5.11) 


The  heat  balance  on  the  bubble- 1 iqu id  interface  reduces  at  this 
order  of  approximation  to  the  following  normalized  equation 

2 


3T 


p  L  Lb 
v  o 


3r 


K  t  (T  -T,  ) 
o  b 


o 


(5.12) 


5.2.3  Interactions 


a.  Order  e 


In  the  asymptotic  theory  presented  here,  the  local 
pressures  and  temperatures  driving  the  growth  of  any  bubble 
are  a  perturbation  of  the  imposed  far  field  pressure, 

P«(t),  and  temperature,  T,.  Since  these  perturbations  are 
due  to  the  presence  of  the  other  bubbles  in  the  flow  field,  the 
leading  terms  can  be  obtained  directly  once  the  first  order 
behavior  of  all  the  bubbles  in  the  cloud  is  determined.  For 
instance,  once  equations  (5.7),  (5.10),  and  (5.12)  are  solved, 

the  variations  with  time  of  the  radius,  ao^(t)  ,  of  any  cavity 
in  the  cloud  are  determined.  This  allows  the  determination  of 
the  intensity  of  all  sources  q0J(t): 


qi<t)  =  Sj  2  P 


(5.13) 
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Consequently ,  the  resultant  "outer"  potential  flow  is  deter¬ 
mined  to  this  order  by: 


♦0(Mr  t) 


(5.14) 


where  M  is  a  field  point,  and  Bl  the  center  of  the  bubble 
B(j),  and  also  the  location  of  the  source  (j).  The  asymptotic 
expansions  of  $0(M,t),  when  the  normalized  distance 

| MB1 |  =  e  r1  goes  to  zero,  contain  additional  terms  other 
than  the  leading  source  term,  qQVr^'  corresponding  to  the 
order  zero  "inner"  potential  flow, 

*o  =  qo(t)/rl  ‘  (5.15) 


These  terms  express  the  interactions  and  are  responsible  for  the 
flow  and  bubble  shape  corrections.  For  instance,  by  application 
of  the  matching  principle  (n  -  m  rule.  Van  Dyke,  1964),  the 
order  e  term  will  determine  the  boundary  condition  at  infinity 
for  the  order  e  "inner"  velocity  potential,  i.e., 

I 

lim  *i  =  l  (  -?.)  ,  (5.16) 

r'i  +  „  3  *o 

where  is  the  initial  distance  between  the  two  cavities’ 

centers  B*  and  B^ . 


In  addition  to  the  at-infinity  boundary  condition  (5.16) 
the  first  correction,  ♦i1,  of  the  undisturbed  potential  flow, 
Go1,  has  to  satisfy  the  Laplace  equation,  (2.1),  as  well  as 
boundary  conditions  on  the  surface  of  the  bubble  B^1).  These 
are  the  contributions  to  order  e  of  the  expansions  in  powers  of 
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e  of  Conditions  (2.2)  and  (2.3)  made  dimensionless.  Similarly, 
the  first  correction,  Tj ,  of  T0  has  to  satisfy  the  equations 
derived  from  (2.5)  and  (2.7). 

Due  to  Condition  (5.16)  the  dynamic  problem  remains 
spherical.  To  this  order  the  effect  of  the  other  bubbles  does 
not  introduce  any  asymmetries,  and  only  changes  the  level  of  the 
velocity  potential.  Therefore,  the  correction  at  this  order 

stems  from  a  modification  in  the  "inner"  problem  of  the  pressure 
imposed  at  infinity  by  the  time  derivative  of  the  added 
at-infinity  velocity  potential,  (5.16).  As  a  result,  the 
solution  of  the  dynamic  problem  at  order  e  is  given  again  by  a 
source  term  which  corrects  the  leading  term  This 

solution  can  be  written: 


where  the  source  intensity,  qj1,  is  given  by: 


i  i2  *i  .  _  l  l 

q  i  =  a >  +  2a^  a„  a 

1  o  1  o  o  j 


(5.18) 


In  order  to  satisfy  the  boundary  conditions  at  the  bubble 
surface,  the  first  correction,  a*1,  of  the  bubble  radius  has 
to  satisfy  the  following  differential  equation,  where  the 
superscript  i  has  been  omitted: 


Here  wjtt)  is  a  correction  of  *0(t)  and  expresses  the  second 
approximation  of  the  value  of  the  vapor  pressure  at  the  bubble 
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wall.  Using  the  expansions  of  the  temperature  in  powers  of  e  as 
in  (5.6),  n0(t)  and  » 1  ( t )  can  be  expressed  as 

ffo(t)  =  [Pv(To(ao't))  "  Pv(T«^  7  ^  '  (5. 20) 

3T  dp 

»l(t)-[pv(T1(a0,t))  +  a,  3^  (ao,t)  ^  (Tq  (aQ  ,t)  )  ]/AP.  (5.21) 

For  the  study  of  the  heat  problem  it  is  useful  to  introduce 
the  following  variable  (again  omitting  the  superscripts  i): 

y  =  ^  [r 3  -  R3  (  6,t)  ]  ,  (5.22) 


by  analogy  with  the  spherical  bubble  case  (Prosperetti  and 
Plesset,  1978).  With  this  variable  change,  the  normalized 
energy  equation  can  be  written: 


T  +  Mil 
3y  1  3r 


1 


4  3T 


-  Rzk]  =  V  1  <r'  17  >  + 


o  , 

r  sin  0 


1  Te  (sln0  H  >  -  7a  (r2  77  sine  7?  >  ] }  +  °U3>.<5.23) 


30 


30 


30 


3y 


After  replacing  r  by  its  value  derived  from  (5.22),  and 
accounting  for  the  expansions  of  R(0,t),  we  obtain  at  the  orders 
e*  and  e,  the  relations: 


and 


0 


(5.24) 


Ti 

where 


-  1 
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e 


3_ 

3y 


(  n 


3T  j 
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3_ 
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3T 
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,(5.25) 
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n  =  (a^  +  3y ) 1/3 


(5.26  ) 


Similarly,  the  heat  balance  condition  on  the  bubble  wall 
becomes  at  order  e°  and  e: 


3T 

»2  — £ 
*o  3y 


y=0 


o  ' 


(5.27) 


and 


*■  3T  i  a  j  3T 

ao  '  3y  ~  3y  n 

1  ao  y=0 


)  -jtzl 


(5.28) 


b.  Higher  Orders 


Continuing  the  same  procedure  as  in  the  preceding  section, 
one  can  derive  the  successive  equations  for  the  flow  field,  the 
temperature  field  and  the  bubble  motion.  The  solution  of  the 
problem  is  made  easier  by  the  use  of  series  expansions  of  the 
velocity  potential  in  spherical  harmonics,  and  of  the  bubble 
surface  equation  and  the  temperature  field  in  Legendre 
polynomials,  Pn(cos0).  The  boundary  condition  at  infinity  for 
any  particular  "inner"  problem  (i),  obtained  by  expanding  the 
expression  of  $  near  B*  (Chahine  and  Bovis,  1983;  Chahine, 

1982)  can  then  be  shown  to  be  up  to  order  e^: 


lim  *i(M,t)  =  l  f-e  X..q^  -  e2  (  cos  e13  +  X.  •  q^ )  + 

i*j  J  J  J 

Jf  >  00  J 

-e3[X?.q3r2  P  2(cos  0i3  )+  X2.  q^  r  cos6i3  +  X..  q3  ]  +  ...}.  (5.29) 
13  1  2 


& 

/v 


Here  qn3  is  the  correction  at  order  en  of  q 03 ,  the 
strength  of  the  source  representing  the  first-approximation 
spherical  volume  change  of  the  bubble  Bw),  and 
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*ij  ■<*</#>  •  (5-30) 

The  superscript  j  denotes  quantities  corresponding  to  the  other 

bubbles,  b(3),  and  9*3  is  the  angle  MB*B3  between 

B*,b3,  and  a  field  point  in  the  fluid,  M,  (see  Figure  1). 

Expressed  in  physical  terms  (velocities,  pressures),  the 
boundary  condition  (5.29)  indicates  that  the  order  e  correction 
to  the  nonperturbed  spherical  flow  field  around  the  test  bubble 
is  a  spherical  modification  of  the  collapse  driving  pressure. 
This  introduces,  as  we  have  seen  in  the  preceding  section,  a 
spherical  correction  ai*(t)  to  the  radius  variations 
a0*(t).  At  the  following  order,  e2,  a  second  correction  of 
the  at-infinity  uniform  pressure  appears,  and  a  uniform  velocity 
field  expressing  a  slip  velocity  between  the  bubble  and  the 
surrounding  fluid  is  to  be  added.  Going  through  the  expansions 
of  the  boundary  conditions  at  the  bubble  surface,  one  can  show 
that  this  induces  a  spherical  correction,  a2*(t),  of  a0*(t), 
and  a  nonspherical  correction  f 2*( t)  »cos 0*9  (Chahine  and 
Bovis,  1983,  Chahine,  1982).  0*9  is  an  angle  which  can  be 
compounded  from  all  the  0*3' s  (see  Figure  3.1).  Things  become 
more  complex  at  the  order  of  expansion  e  ,  where  in  addition  to 
the  uniform  pressure  and  velocity  corrections,  a  velocity 
gradient  generated  by  the  flow  field  associated  with  the  motion 
of  all  the  other  bubbles,  is  to  be  accounted  for  to  generate  a 
nonspherical  correction  of  form  P2(cos0*9). 

Resulting  from  the  above  remarks  on  the  at-infinity 
boundary  condition,  one  can  show  that  the  equation  of  the 
surface  of  the  bubble  B  (  * )  and  the  temperature  can  be  written 
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Rl(6ig,t)  =  a*(t)  +  e  a  i  (t)  +  e2  [a2(t)  +  f2  (t)  •  cos0ig]  + 

+  e3  [a  3<  t)  +  f  3  ( t )  *cos  0 ig  +  gj(t)  «P  2  (cos9lg)]  +  o(  e3)  , 

(5.31) 

T(r, 0ig ,t)=TQ (r,t)  +  e  Tj(r,t)  +  e2  [T 2 0{ r , t ) +T 2  i ( r , 1 )  *cos 9ig  ]+ 

+  e  3  [T  3  o  ( r ,  t )  +  T31(r,t)  •cos9ig+T32(r,t)  •  P  2(cos  9ig )  ]+o(  e  3)  , 

(5.32) 


provided  that  the  initial  bubble  shape  is  spherical.  Therefore, 
up  to  the  order  e3  each  inner  problem  is  axisymmetr ic ,  and  the 
axis  of  symmetry  for  every  bubble  is  in  the  direction,  B^G,  of 
its  motion  towards  the  bubble  cloud  "center". 

We  introduce  dn  defined  as  the  sum  of  the  deformation 

rate,  f n ,  and  of  the  origin  of  axis  translation  velocity, 

» 

bn,  both  of  form  cos  9: 


d  =  £  +  b'  . 

n  n  n 


(5.33) 


One  then  obtains  the  following  differential  equations  for  the 
order  e2  radius  components,  and  d2,  similar  to  (3.11). 


••  2  .  •  • 

aoa2  +  3  50^2  +  a2  (ao  "  2We_  1  ao~2) 


•  -  —  i  _  —  3 

-  %  a 2  -  a  i  a  j  -  2w  a2  a 
2  i  11  e  i  o 


-~3  •  . 

2  a  +  n2(t)  -  l  X  q\  ,  (5. 
1  i*i  J 


a0d2  +  3  aj,  =  -  3  I  q  <&„<£  +  a'o43) 


cos  9 


(5.35) 


wV  v  • "*'■  - 


.‘s’*  /•  .’* ' v  ■/  *.*  v 
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In  the  above  equations,  the  superscript  i  has  been  omitted, 
and  ao<t)  and  a^t)  are  obtained  at  the  preceding  orders  by 
solving  the  differential  equations  (5.7)  and  (5.19).  Hnm(t) 
are  higher  order  approximations  of  the  nondimensional  difference 
between  the  vapor  pressure  at  any  time  and  its  initial  value. 

The  indexes  n,m  correspond  to  those  used  in  the  expansions  of 
the  temperature.  Equation  (5.32). 

The  energy  equations  at  order  e2  can  be  written  using 


ft  =  (  a  o  3  +  3y) 


l/  3 


(5.36) 


i  1  3T 

T  2  o  -  U4  -ay2  ) 


3T 


3T  ■ 


>  ~r—  r  ( 4  a  a  2n  +  4a2a?n  +  2a4a2n  2)  -r-°  +  4a2ain  — — 

e  3y  LV  o  i  o  2  o  i  3y  o  1  3y 


(5.37) 


/V/ 

T  -  P 
2  1  e 


i  a  3T  2 1 

3y  (  n  3y 


2h 


~27 


)  -  2rf2  T21  ]  = 


cos9  J  3T 

( -  +  2&  a  f2  +  a2f  2  +  l  X2 .  q-1  - . —  n2)  + 

n  o  o  *  o  *  .  lj  ^o  ,ig  ' 

J  cos  e  ^ 


3y 


-1  r?  ~2  I2  f  3To  .  3  ,,  2  ”2  3To 

l2ri  ao  f2  IT-  +  177  (4r  f  2  ar 


)]  .  (5.38) 


e  '  ~o  3y  3y  '  -  "  -  i  -0  3y 

The  corresponding  conditions  of  heat  balance  on  the  bubble  wall 
are 


t  ay 


3T  20  2a  j  3T  j 

+  -  +  2  ( 


a  o 


3y 


-2 

a  2  a  j  3T 

+  —  >  ly2  Jy.0****4*  '(5-39) 


a  o  a 
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3T21  f2 

[  - r - +  2  -  + 

a0 


3T 
_ O 

3y 


]y=0=^f2 


(5.40) 


Similar  equations  can  be  derived  for  order  e3  (Chahine  and  Liu, 
1983)  and  are  not  presented  here  for  conciseness. 

5.3  Numerical  Solution 


5.3.1  General  Solution 


The  system  of  equations  derived  up  to  order  e3  constitutes 
a  set  of  14  equations  for  the  14  unknown  components  of 
Ri(6,t),  and  T*(0,t)  (expansions  (5.31)  and  (5.32)).  By 
solving  this  system  one  determines  completely  the  flow  and 
temperature  fields  as  well  as  the  bubble  motion  and 

deformation.  A  numerical  solution  of  these  equations  is 
feasible  and  could  be  performed  using  the  same  procedure  as 
Dalle  Donne  and  Ferranti  (1975).  Their  study  dealt  with  a 
single  bubble  growth  and  thus  solved  only  equations  (5.7), 
(5.10),  and  (5.12).  Here  the  same  approach  would  have  to  be 
performed  for  all  seven  components  of  the  bubble  radius  (up  to 

e3). 

Since  the  equations  are  not  independent,  the  procedure 
would  start  by  determining  at  a  given  time  step  the 
temperature  at  the  lowest  order  of  approximation  (e°)  and  the 
corresponding  radius  approximation.  Knowing  this,  one  can 
compute  the  successive  temperature  corrections,  and  the 
successive  radius  corrections.  At  each  time  step  an  iteration 
process  would  be  used  to  insure  a  good  correspondence  between 
the  obtained  temperature  and  radius  values.  Stepping  in  time  of 
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the  computation  could  be  obtained  with  a  Runge-Kutta  procedure 
which  solves  each  differential  equation  yielding  the  bubble 
radius  values.  The  determination  of  the  temperature  field  is 
more  elaborate  and  requires  a  stepping  in  both  time  and  in 
space.  This  latter  involves  writing  a  finite  difference  scheme 
and  replacing  the  integration  field  with  a  grid  of  mesh  points. 
This  general  solution  is  not  developed  here.  We  considered 
instead  the  cases  where  large  initial  superheats  make  a  thermal 
boundary  layer  approximation  valid. 

5.3.2  Thermal  Boundary  Layer  Approximation 

If  the  distance  <5  in  which  the  temperature  rises  from  its 
value  at  the  bubble  wall  to  approximately  the  imposed  ambient 
temperature,  T„,  is  small  compared  with  the  bubble  radius,  R, 
an  approximate  solution  can  be  obtained  more  easily  than  with 
the  method  described  in  the  preceding  paragraph.  By  considering 
heat  diffusion  in  the  liquid,  spherical  bubble  growth  rate,  and 
a  heat  balance  at  the  bubble-liquid  interface,  Plesset  and 
Prosperetti  (1977)  estimate  S/R  by: 

,  D  L  p  p  y 

_®  ..  _ Y  —  _ Y  _ h _  _  1  /  r  41% 

R  K  (Tb  -  Tb)  '  p  C(T#-  Tb)  *  U 

Thus,  a  boundary  layer  approximation  is  valid  as  long  as  the 
Jacob  number,  J,  is  much  larger  than  one.  For  a  spherical 
bubble,  comparisons  between  numerical  computations  obtained 
using  this  approximation  and  those  obtained  by  solving  the  exact 
equations  gave  very  close  agreement  for  J  >_  3,  (Prosperetti  and 
Plesset  1978  and  Plesset.  1980). 
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When  the  boundary  layer  approximation  is  used,  the  system 
of  heat  equations,  presented  above,  simplifies  considerably. 
Indeed,  in  that  case  the  temperature  departs  from  T„  only  in 
the  liquid  region  close  to  the  bubble-liquid  interface,  and  the 
values  of  r  which  are  of  interest  are  close  to  R(0,t).  The 
variable,  y,  defined  in  (5.22),  is  then  small  compared  to  a3, 
and  we  can  write 

y  =  £  Y  aQ  3  »  (5.42) 

where  y  and  aQ  are  of  order  1,  and  £  is  a  small  parameter 
[  5=0 ( J“  1 )  ]  •  The  problem  considered  then  contains  two  small 
parameters  e  and  5,  and  an  asymptotic  solution  uniformly  valid 
when  both  e  and  £  go  to  zero  can  be  obtained  when  a 
relationship  between  the  two  parameters  is  defined  through  the 
use  of  the  principle  of  least  degeneracy  (Darrozes,  1971). 

Considering  the  heat  equation,  (5.23),  one  can  determine 
the  relation  needed  between  e  and  £  to  conserve  the  maximum 
number  of  terms  in  the  leading  orders  of  approximation.  In 
order  to  prevent  the  order  z  0  expansion,  (5.24),  from 
degenerating  when  5  goes  to  zero,  the  Peclet  number  has  to  be 
large  enough  to  satisfy 

Pe  =  0(  T2)  ,  (5.43) 

in  which  case  both  terms  of  the  equation  are  conserved. 
Similarly,  to  conserve  the  maximum  terms  at  the  following  order, 
e,  one  needs  to  keep  together  the  leading  terms  coming  from  the 
expansions  in  powers  of  e  and  those  from  expansions  in  power  of 
£  (e.g.,  in  the  expansions  of  n 4 )  -  This  "least  degeneracy"  is 
obtained  when 


e  =  o(  e) 


(5.44) 
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Using  the  relationships  (5.43)  and  (5.44)  between  Pe ,  C,  and 
e,  the  expansions  become  straightforward. 


The  equations  obtained  at  the  first  order  expansion  in  both 
parameters  (orders  e#  and  5°)  are  those  for  the  case  of  an 
isolated  bubble.  These  equations  are  the  Rayleigh-Plesset 
equation  (5.7)  and  the  heat  equation  (5.24)  in  which  n4  takes 
the  value  a,;,  and  which  becomes: 


+  p 


-  l 


(5.45) 


A  solution  is  readily  available  for  this  case  and  was  derived  by 
Plesset  and  Zwick  (1952)  and  Forster  and  Zuber  (1954)  using 
Laplace  transform  methods.  The  nondimensional  temperature  at 
the  bubble  wall  is  given  by: 


[  /fca 4(y)  dy] 
x  u  (5.46) 


where,  to  be  consistent  with  the  assumptions  made  in  deriving 
this  solution,  D  and  K  are  constant  and  evaluated  at  T5  while 
L  and  pv  are  functions  of  time.  The  numerical  procedure  is 
greatly  simplified  now  that  an  analytical  expression  for  the 
temperature  at  the  bubble  wall  is  known.  The  finite  difference 
method  which  would  have  been  used  in  the  general  case  is  here 
replaced  by  a  numerical  computation  of  the  integral  equation 
(5.46).  An  iteration  procedure  is  required  to  insure  that  the 
computed  value  of  T0(a0,t)  does  not  differ  significantly 
from  the  value  presumed  in  the  computation  of  the  integrand. 


Tracor  Hydronautics 


-73- 


Plesset  and  Zwick  (1952)  also  gave  the  solution  of  the 
problem  when  equation  (5.46)  contains  a  right  hand  side  which  is 
a  known  function  of  time  (heat  source  term).  Using  a  matched 
asymptotic  procedure  they  also  computed  the  following  order  of 
approximation,  0(5).  These  solutions  correspond  to  the 
following  order  equations  in  powers  of  £  and  5  for  the 
multibubble  problem.  For  simplicity  these  equations  will  not  be 
listed  here. 

5.4  Numerical  Examples  for  Symmetrical  Clouds 

In  order  to  illustrate  the  method  presented  above  we 
consider  numerical  solutions  for  a  cloud  of  simple  geometry. 

The  bubbles  are  distributed  in  a  symmetrical  configuration  and 
are  initially  of  equal  size.  With  this  configuration  all 
bubbles  have  the  same  radius  history.  All  summations  in  the 
dynamic  equations  (5.19,  5.34,  5.35)  reduce  to  multiplications 
of  the  characteristics  of  a  single  bubble  by  one  of  the  three 
geometrical  constants  clf  c2*  c3  defined  earlier  in  (3.13). 

An  additional  simplification  of  the  numerical  solution  can 
be  introduced  if  one  notices  that  during  the  bubble  growth  the 
departure  from  the  initial  spherical  shape  happens  very  late  in 
the  bubble  history  and  only  when  the  asymptotic  approach  starts 
losing  its  validity.  This  is  not  true  for  the  cloud  collapse 
(see  preceding  section).  Figure  5.1  shows  the  variation  with 
time  of  the  major  radius  of  an  individual  bubble  in  a  cloud 
configuration  of  N  bubbles  symmetrically  located  on  a  sphere. 

For  this  figure,  heat  transfer  has  been  neglected.  We  observe, 
for  the  isolated  bubble,  the  well  known  asymptotic  linear  growth 
behavior.  However,  when  the  number  of  interacting  bubbles 
increases,  the  pressure  field  associated  with  the  dynamics  of 
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the  other  bubbles  in  the  cloud  reduces  the  growth 
rate  of  the  test  bubble.  This  deviation  increases  with  the 
number  of  bubbles,  N,  until  for  N  =  12  for  the  case  studied,  the 
method  apparently  fails  for  t  >  0.1.  The  radius  corrections 
(illustrated  in  the  figure  by  the  amount  of  deviation  of  the 
radius  in  an  N-bubble  case  from  the  isolated  bubble  case)  become 
large  compared  with  the  order  zero  radius.  Figure  5.2  shows, 
for  the  same  bubble  configuration,  the  ratio  of  the  nonspherical 
to  the  spherical  part  of  R(0,t)  in  the  expansion  (5.31).  In  all 
cases  but  the  obvious  one  where  the  method  breaks  down,  the 
relative  deformations  remain  less  than  4  percent  while  the 
bubble  radius  is  2,000  times  its  initial  value.  Based  on  this 
observation  and  as  a  first  step  towards  a  more  precise  solution, 
we  have  neglected  in  the  numerical  program  developed  the 
contribution  of  nonsphericity  to  the  heat  transfer  problem. 
Therefore,  the  temperature  field  was  approximated  by  a 
spherically  symmetrical  field.  However,  this  field  accounts  for 
interactions  and  differs  from  that  of  the  isolated  bubble  case 
because  of  the  contributions  of  the  higher  order  spherical  terms 
of  the  bubble  equation.  Indeed,  Equation  (5.46),  relating  the 
bubble  wall  temperature  to  a  spherical  bubble  radius  history, 
was  applied  to  the  spherical  part  of  the  bubble  radius,  i.e., 
to 

A(t)=a0(t)  +  e  3l(t)  +  e2  a2(t)  +  e3  a3(t)  .  (5.47) 

With  this  simplification,  at  any  time  step  all  dynamical 
equations  are  solved  using  the  value  of  the  vapor  pressure 
corresponding  to  the  liquid  temperature  at  the  radial  distance 
A ( t ) .  This  temperature  is  computed  at  the  preceding  time  step 
using  equation  (5.46).  The  nonspherical  part  of  the  bubble 
shape  is  not  disregarded  and  is  computed  neglecting  variations 


Tracor  Hydronautics 


of  the  liquid  temperature  along  the  bubble  surface.  This  is 
valid  as  long  as  the  bubble  deformation  is  negligible.  Since 
we  restrict  this  study  to  that  case,  the  validity  of  the  results 
is  checked  by  monitoring  the  relative  value  of  the  computed 
nonspherical  to  the  spherical  components  of  the  bubble  surface 
equation.  The  computation  is  stopped  when  an  imposed  limit  is 
exceeded.  Figure  5.3  presents  a  flow  chart  of  the  developed 
code  HOTCLD. 

A  series  of  numerical  cases  was  studied  using  a  VAX  11/750 
computer.  We  have  considered  different  variations  of  the  number 
of  bubbles  and  configuration,  the  ambient  pressures,  the  initial 
bubble  radius,  and  the  amount  of  superheat.  The  duration  of  a 
typical  run  was  about  10  minutes  of  CPU  time  (for  2,000  time 
steps).  The  computation  involves  the  solution  of  the  heat  and 
dynamical  equations  for  an  N  bubble  configuration,  the  study  of 
the  corresponding  case  of  an  isolated  bubble  with  and  without 
heat  transfer,  and  the  computation  of  pressure  histories  at 
three  locations  in  the  flow  field. 

In  all  the  presented  figures,  the  curves  are  stopped  when 
the  computations  become  invalid  due  to  large  bubble 
interactions.  Figure  5.4  shows  clearly  the  influence  of 
interactions  on  the  bubble  radius  history.  Since  the  bubble 
does  not  remain  spherical,  the  value  of  Rie1^)  represented  in 
this  figure  corresponds  to  the  point  on  the  bubble  closest  to 
the  cloud  center,  the  "lower-minor  radius.”  (See  Figure  3.1.) 
The  classical  results  of  asymptotic  growth  in  t  for  the 
inertia-controlled  bubble  expansion  and  in  t°  for  the 
heat-controlled  bubble  expansion  can  be  seen.  If  there  was 
no  pressure  drop  a  would  be  1/2.  However,  here  a  is  much  closer 
to  1 ,  as  obtained  by  earlier  studies  on  single  bubbles  (Jones 
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and  Zuber,  1978;  Theofanous  et  al . ,  1969;  and  Cha  and  Henry, 
1981).  The  most  important  result  obtained  here  is  that  bubble 
growth  is  inhibited  by  bubble  interactions.  Very  clearly  at  a 
given  time  the  bubble  size  decreases  with  the  number  of 
interacting  bubbles.  This  decrease  exceeds  20  percent  for  a  5 
bubble  system  for  nondimensional  times  larger  than  10,  or  one 
millisecond  after  the  start  of  the  growth  (case  of  Figure  5.4). 

Figure  5.5  shows  the  effect  of  bubble  interactions  on  the 
liquid  temperature  at  the  bubble  wall.  The  presence  of  other 
growing  bubbles  in  the  field  is  seen  to  reduce  the  heat  trans¬ 
fer  at  the  bubble  wall  and  thus  the  temperature  drop  in  its 
vicinity.  For  example,  for  a  five-bubble  system  the  deviation 
from  the  isolated  bubble  case  of  the  temperature  drop  is  more 
than  30  degrees  one  millisecond  after  the  initial  pressure 
drop.  This  result,  coupled  with  that  obtained  for  the 
variations  of  the  bubble  radius,  is  important  for  any  practical 
computation  of  heat  transfer  in  a  two-phase  medium. 

Figure  5.6  shows  the  modification  of  the  bubble  shape 
during  its  growth  for  the  same  N-bubble  systems  shown  in  Figures 
5.4  and  5.5.  Represented  are  the  bubble  shapes  at  two  instants 
during  the  growth  process.  As  expected,  in  the  presence  of  an 
N-bubble  cloud,  the  side  of  the  bubble  facing  the  cloud  center 
is  seen  to  be  slightly  "pushed  away"  from  the  cloud  center  and 
the  bubble  is  seen  to  elongate  in  a  direction  tangential  to  the 
sphere.  However,  any  point  on  its  surface  always  remains  inside 
the  corresponding  fictitious  isolated  bubble  growing  under  the 
same  conditions.  The  deformation  decreases  as  the  number  of 
interacting  bubbles  increases. 
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Figures  5.7  and  5.8  show  the  influence  of  the  amount  of 
pressure  drop  and  initial  superheat  on  the  bubble  growth.  These 
figures  consider  an  isolated  bubble  as  well  as  a  five-bubble 
system.  The  same  remarks  made  in  the  preceding  paragraphs  apply 
here  when  the  influence  of  the  number  of  bubbles  is  considered. 
In  all  cases  the  initial  bubble  radius  is  the  same,  and  the 
pressure  drops  to  the  same  value.  However,  since  the  initial 
pressures  vary  from  one  case  to  another  and  since  all  bubbles 
are  considered  to  be  initially  at  equilibrium,  the  initial 
temperature  and  thus  the  initial  amount  of  superheat  varies  from 
one  case  to  another.  To  isolate  the  two  effects  one  has  to 
consider  the  case  where  the  bubbles  are  not  initially  at 
equilibrium.  Another  option  would  be  to  have  the  same  initial 
pressure,  radius  and  temperature  and  to  vary  the  value  of 
pinf •  We  consider  this  case  below.  Figures  5.7  and  5.8  show 
that  the  normalized  bubble  radii  and  growth  rates  are  larger  at 
any  given  time  when  the  amounts  of  pressure  drop  and  superheat 
are  greater.  In  the  absence  of  heat  transfer,  scaling  effects 
are  mainly  due  to  the  differences  in  the  Weber  number,  We  and 
the  initial  pressure  parameter  .  When  heat  transfer  effects 
are  included,  there  is  an  additional  parameter,  the  Jacob 
number.  These  effects  counterbalance  each  other  in  real  time, 
and  one  observes  a  minor  influence  of  the  initial  value  of  the 
pressure  (for  the  same  initial  radius)  when  the  radius 
variations  are  plotted  with  dimensional  variables  (Chahine  and 
Liu ,  1983). 

A  similar  result  is  seen  when,  for  the  same  initial  bubble 
radius  and  liquid  temperature,  the  ambient  pressure  drops  from 
the  same  initial  pressure  to  different  subsequent  values.  In 
this  case  the  initial  amount  of  superheat  is  the  same  for  all 
the  cases  of  pressure  drop  studied.  However,  the  subsequent 
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amounts  of  superheat  differ  from  one  case  to  another.  The  use 
of  nondimensional  variables  decreases  (but  does  not  cancel, 
because  of  nonlinearities)  the  influence  of  dynamic  factors, 
however  it  does  not  alter  thermal  effects.  The  results  (Figure 
5.9)  show  again  a  larger  bubble  radius  during  the  growth  period 
for  higher  pressure  drops.  The  comparison  with  the  five-bubble 
case  can  also  be  observed  and  shows  again  the  retarding  effect 
on  bubble  growth  due  to  collective  bubble  behavior.  Figure  5.10 
completes  the  picture  by  showing  the  temperature  drop  at  the 
bubble  wall  for  the  different  cases  studied.  The  same 
observations  as  those  made  above  are  repeated,  namely  higher 
temperature  drops  for  smaller  pressure  drops  or  higher  number  of 
interacting  bubbles. 

The  last  series  of  results  consider  the  influence  of  the 
initial  bubble  size  for  given  fixed  pressure  conditions.  With 
the  assumption  that  the  bubble  is  initially  at  equilibrium,  the 
modification  of  the  initial  bubble  size  also  corresponds  to  a 
change  of  the  amount  of  superheat.  Figure  5.11  shows  the  pre¬ 
dominance  of  the  effect  of  the  amount  of  superheat  factor  on  the 
bubble  growth;  initially  smaller  bubbles  attain  greater  sizes 
because  of  larger  amounts  of  superheat.  This  effect  is, 

however,  coupled  with  the  nonlinearities  of  the  dynamical 
equations  which  favor  smaller  initial  bubble  radii  in  the  first 
phase  of  the  growth.  Figure  5.12  shows  the  same  effect  with 
nondimensional  variables  and  compares  a  five-bubble  system  with 
the  isolated  bubble  case.  One  can  notice  that  the  inhibition 
effect  due  to  bubble  interactions  is  larger  for  smaller  initial 
bubbles  or  larger  amounts  of  superheat.  Finally,  Figure  5.13 
describes  the  temperature  drop  at  the  bubble  wall  for  the  same 


cases. 


i 
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The  results  obtained  show  that  a  significant  influence  of 
bubble  interactions  on  bubble  growth  and  heat  transfer  exists. 
The  effects  of  this  influence  can  be  summarized  as  follows: 

a.  The  growth  rate  of  the  bubbles  is  reduced, 

b.  The  radius  of  any  bubble  at  a  given  time  is  smaller 
than  would  be  found  for  an  isolated  bubble,  and 


c.  The  temperature  drop  at  the  bubble  wall  is  smaller  at 
any  given  time  than  would  be  found  for  an  isolated  bubble. 

These  effects  increase  with  the  number  of  interacting 
bubbles  as  well  as  with  the  amount  of  superheat  and  pressure 
drop.  These  results,  which  were  obtained  using  small 
perturbations  assumptions,  are  expected  to  remain  valid  and 
become  more  significant  when  the  void  fraction  becomes  larger. 
Accounting  for  these  effects  is  important  for  increasing  the 
accuracies  of  the  existing  transient  two-phase  flow  codes. 

The  study  presented  here  could  be  improved  by  introducing  a 
finite  speed  wave  propagation  in  the  cloud  and  by  accounting  for 
the  compressibility  of  the  medium.  The  analytical  equations 
derived  for  the  general  bubble  configuration  case  (no  symmetry 
and  unequal  bubble  size)  could  be  expanded  to  a  numerical 
approach  in  a  relatively  simple  manner.  The  resolution  of  the 
problem  could  also  be  extended  to  low  superheat  cases  (small 
Jacob  numbers)  and  to  larger  interactions  and  bubble 
deformations  by  numerically  implementing  the  analytical  approach 
presented  above  which  was  not  used  in  the  numerical  examples. 
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6.0  INCLUSION  OF  GAS  DIFFUSION  AT  THE  BUBBLE  WALL 

6.1  Introduction 


An  analysis  was  developed  for  the  inclusion  of  the  effects 
of  gas  diffusion  into  and  out  of  cavitation  bubbles  in  clouds. 
The  solution  scheme  is  similar  to  that  used  by  Chahine  and  Liu 
(1985)  for  heat  transfer  to  a  bubble  in  a  superheated  liquid 
(see  Section  5).  The  equations  governing  mass  diffusion  have 
been  coupled  to  those  governing  bubble  dynamics  and 
non-spher ical  deformations.  The  equations  solved  include  the 
Rayleigh-Plesset  equation  for  bubble  dynamics,  the  transport 
eqration  for  the  dissolved  gas  in  the  liquid,  and  an  energy 
equation  for  the  bubble  (first  law  of  thermodynamics).  An  ideal 
gas  equation  of  state  is  assumed  for  the  bubble  which  is  taken 
to  contain  a  mixture  of  noncondensible  gas  and  vapor.  These  two 
components  are  considered  to  form  an  ideal  gas  mixture.  Due  to 
the  relatively  very  short  characteristic  time  of  vaporization 
compared  to  bubble  dynamic  and  gas  diffusion  characteristic 
times,  the  quantity  of  vapor  in  the  bubble  is  assumed  to  vary 
instantaneously  such  that  the  vapor  pressure  remains  constant 
and  at  the  equilibrium  value  for  the  liquid  temperature.  A 
boundary  layer  solution  analogous  to  that  of  Plesset  and  Zwick 
(1952)  is  utilized  for  the  gas  transport  equation.  This 
provides  an  expression  for  the  gas  concentration  at  the  bubble 
surface  as  a  function  of  time.  The  concentration  of  dissolved 
gas  at  the  bubble  surface  is  related  to  the  partial  pressure  of 
the  gas  within  the  bubble  by  Henry's  law.  The  gas  concentration 
varies  from  the  ambient  initial  value  in  the  liquid  to  the 
bubble  wall  value  over  a  thin  boundary  layer. 
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6 . 2  Equations  Governing  Mass  Diffusion 
6.2.1  General  Equations 

The  transport  of  noncondensible  gas  in  the  liquid  medium  is 
governed  by 

H  +  *  •  VC  =  Dg  J2C  (6.1) 

subject  to  the  boundary  and  initial  conditions: 

C  -*•  C  for  r  ®  , 

00 

C  =  Cb  for  t  0  and  r  R  ,  (6.2) 

C  «  Csat '  r  "  R'  >  0  ‘ 

Here,  C  is  the  concentration  of  dissolved  gas  in  the  liquid  in 

units  of  moles  per  unit  volume,  Dg  is  the  molar  diffusivity  of 

the  gas  component  in  the  liquid,  and  R  is  the  radial  location  of 
the  bubble  wall.  In  general,  R  =  R(  9,  <J>,  t).  The  gas  concentra¬ 
tion  is  taken  to  be  saturated  at  the  bubble  wall:  C  =  Csaf 
With  this  assumption,  the  gas  concentration  at  the  bubble  wall 
can  be  related  to  the  partial  pressure  of  that  gas  component  in 
the  bubble  by  Henry's  Law: 

pg  =  H  Csat  =  H  •  C( r=R)  ,  (6.3) 

where  Pg  is  the  partial  pressure  of  noncondensible  gas  within 
the  bubble  (taken  to  be  uniform  throughout  the  bubble)  and  H  is 
Henry's  constant  at  the  liquid  temperature,  a  property  of  the 
gas-liquid  combination. 
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As  the  bubble  grows  or  shrinks,  noncondensible  gas  will  be 
transported  across  the  bubble  wall.  The  net  transport  of  gas 
into  the  bubble  can  be  related  to  the  gradient  of  gas  concentra 
tion  in  the  liquid  at  the  bubble  wall: 


D 


g 


/ 

s 


3C 

3n 


r=R 


dS 


(6.4) 


where  S  is  the  surface  of  the  bubble,  n  is  the  direction  normal 
to  the  bubble  surface,  and  ng  is  the  time  rate  of  change  of 
the  total  number  of  moles  of  gas,  ng ,  within  the  bubble. 


6.2.2  Boundary  Layer  Solution  for  Gas  Concentration 


A  general,  detailed  solution  of  equations  (6. 1-6.3)  would 
likely  involve  a  time  consuming  numerical  approach,  such  as  a 
space  and  time  dependent  finite  difference  scheme,  and  is  not 
necessary  for  many  cases  of  interest.  Plesset  and  Zwick  (1952) 
and  Forster  and  Zuber  (1954)  obtained  a  solution  to  (6. 1-6.3) 
for  the  case  of  an  isolated  spherical  bubble  for  which 
appreciable  concentration  gradients  are  confined  to  a  boundary 
layer  of  thickness  6  which  is  small  compared  to  the  bubble 
radius : 


C  =  C  - 
00 


R2  i£ 

3r 

r  =  R 

[  /*  R4  (Y)  <3y  ]1/2 


dx. 


(6.5) 


This  expression  applies  to  the  first  order  approximation  of  the 
bubble  radius,  R  =  ao(t).  A  similar  relation  is  found 
applicable  to  R  =  aQ(t)  in  the  heat  transfer  study  (see 
Section  5.0  and  Equation  5.46).  Equation  (6.4)  can  then  be 
expressed  as 
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3C 

3r 


n 


r=a. 


4  ir  D 


(6.6) 


and  substituted  into  (6.5)  to  obtain 


C  =  C  -  [  it  d  ]  1/2  fj* 


3 

3x  ng 


•  °  J°  t  J*  a*  (y)  dy'F2' 


dx .  (6.7) 


We  now  nondimensionalize  the  variables  in  (6.7)  and  define 


C  =  C/C 


n  =  n  /n  , 

g  g  gQ 


t,  xf  y  *  t/t,  x/x,  y/T  , 


(6.8) 


5o  =  ao/rb  ' 


r  =  r/rb  , 
o 


where 


C«,  is  the  gas  concentration  for  r  -*■«, 

rbQ  is  a  reference  characteristic  bubble  radius, 

ngQ  is  the  number  of  moles  of  gas  in  a  bubble  of  radius 
rbQ  at  equilibrium  at  t=0,  and 

t  =  rboyp/ AP 1  is  the  characteristic  time  for  growth  of  a 
bubble  of  radius  rb  . 
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Equation  (6.7)  becomes,  in  nondimensional  form: 


C  =  1  - 


[ — - ] 

Ut  Dg 


1/2 


n 


4  irr, 


/  o 


/j  "ao  (*)d^ 


T72 


dx  .(6.9) 


Let  pgQ  be  the  initial  density  of  the  gas  in  the  bubble. 
Then 


n 


m  /  M 
<3  q 

3o  * 


ir  r, 


p  /  M 

V  9 


(6.10) 


where  Mg  is  the  gas  molecular  weight  and  mgo  is  the  mass 
of  gas  in  the  reference  bubble  of  radius  rbQ.  Equation 
(6.9)  can  then  be  written  as 


C  =  1  -  T  [ 


1/2 


ITT  D 


M  C 

g  ■ 


•  r\ 


n 


[  f*  ( y )  dy  ] 

x 


172 


dx,  (6.11) 


for  C„  expressed  in  moles  per  unit  volume,  and  pgQ  in  mass 
per  unit  volume.  If  the  dissolved  gas  concentration  is 
expressed  as  a  volume  fraction  (e.g.,  cm 3  of  gas  per  cm 3  of 
water),  we  will  let  C,*  be  the  concentration  at  infinity  and 

CV  =  C  M  /  d  .  (6 

”  *  g  gQ 


Then  equation  (6.11)  becomes 
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1/2 


CV 


=  1  - 


3C 


IT  T  D 


[  I-  (y)dy] 

X 


T72 


dx  ,  (6.13) 


where  Cv  is  now  a  volume  fraction  concentration 

V 

nondimensionalized  on  C*.  We  will  use  expression  (6.13)  for 
the  transport  of  noncondensible  gas  across  the  bubble  wall. 


6..Z.3  Equation  of  State  Within  the  Bubble 


The  gas  and  vapor  within  the  bubble  are  taken  to  comprise 
an  ideal  gas  mixture,  and  the  equation  of  state  for  that  mixture 
is  given  by 


+  V 


nv> 


(6.14) 


where 


Pg ,  Pv  =  partial  pressures  of  gas  and  vapor, 

Vfc  =  volume  of  bubble, 

ng ,  nv  =  number  of  moles  of  gas  and  vapor  within  the 
bubble , 

Ru  =  Universal  gas  constant,  and 

Tg  =  Absolute  temperature  of  the  gas  and  vapor  mixture 
within  the  bubble. 

The  temperatures  of  the  gas  and  vapor  phases  are  taken  to  be  the 
same,  Tg .  For  an  ideal  gas  mixture,  each  component  also  obeys 
the  ideal  gas  law: 
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P  V.  =  n  R  T  , 
g  b  g  u  g 


Pv  Vb  =  nv  Ru  Tg 


(6.15) 


In  this  analysis,  the  partial  pressure  of  the  vapor  is  fixed  (at 
the  equilibrium  value  for  the  constant  liquid  temperature)  while 
Pg  is  allowed  to  vary  with  bubble  behavior,  mass  transfer,  and 
the  bubble  energy  balance. 


We  nondimensionalize  with  the  characteristic  value  of  pressure 
change,  AP,  and  the  equilibrium  values  at  the  characteristics 
bubble  radius,  rb0: 

P  ,  P  =  P  /AP,  P  /  AP  , 
g  v  g  v 


n  ,  n  =  n  /n  ,  n  /n 
g  v  g  gQ'  v'  gQ 


V, 


-  V<  i  *  O  -  5o 

o 


(6.16) 


The  temperature  is  normalized  on  the  liquid  temperature, 

Tg  =  Tg/T x,  (6.17) 

Equation  (6.14)  can  then  be  written  as 


(Pg  +  Pv>  Vb  =  S  (ng  +  nv»  Tg  '  '6-,8> 


where 
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since 


P 


(6.20) 


the  partial  pressure  of  gas  in  a  reference  bubble  of  radius 
rt>Q  at  equilibrium  for  the  initial  liquid  pressure  and 
temperature.  Thus,  the  quantity  8  can  be  viewed  both  as  a  ratio 
of  the  reference  bubble  gas  partial  pressure  to  driving  liquid 
pressure  change  or  as  a  dimensionless  universal  gas  constant. 
Similarly,  we  write 


6  n. 


T 

g 


(6.21 ) 


Use  of  the  ideal  gas  equations  of  state  has  introduced 
another  unknown,  Tg ,  and  thus  necessitates  use  of  another 
condition.  This  is  provided  by  the  first  law  energy  balance  of 
the  ideal  gas  mixture  within  the  bubble.  Two  cases  are 
considered.  In  the  first  case,  the  approximation  is  made  that 
the  gas  and  vapor  mixture  temperature  is  constant  and  equal  to 
that  of  the  liquid.  This  approximation  is  valid  when  the  bubble 
is  initially  at  the  liquid  temperature,  and  changes  in  bubble 
size  occur  slowly  enough  such  that  heat  transfer  can  equilibrate 
the  bubble  and  liquid  temperatures.  The  second,  more  compli¬ 
cated  case,  involves  writing  the  first  law  energy  balance  for 
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the  bubble  with  mass  transfer  across  its  surface.  This  formula 
tion  allows  for  Tg  to  differ  from  without  an  arbitrary 
adoption  of  a  polytropic  constant.  Its  use  actually  enables 
determination  of  the  effective  polytropic  constant. 

In  the  first  approach,  the  change  in  state  of  the 
isothermal  ideal  gas  is  given  by 


p  v  P  V, 

g  vb  .  g  vb 

—2 -  =  constant  =  — - 


—  =  Pg  ' 

ng  R  -  rb 

o 


=  P 


go  a^ 
o 


(6.22) 


This  provides  a  simple  relation  between  the  number  of  moles  of 
gas,  the  gas  partial  pressure,  and  the  bubble  volume  at 
different  times. 

6.2.4  Energy  Balance 

For  the  second  approach,  we  consider  the  bubble  wall  to 
constitute  a  deformable  and  permeable  control  surface  and  write 
the  first  law  energy  balance  for  the  control  volume  bounded  by 
that  surface: 


d  U  = 


-  d  W  +  J  ni  hi  dt  ;  i  =  v,  g  , 


(6.23) 


where 


=  Change  in  internal  energy  of  the  control  volume  in 
time  dt. 
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dW  =  Work  done  by  the  control  volume  on  its  surroundings  in 
time  dt, 

n^  =  Net  molar  rate  of  mass  transfer  of  component  i  into 
the  control  volume,  and 

hi  =  Specific  enthalpy  of  component  i. 

In  (6.23),  we  have  neglected  changes  in  kinetic  and  potential 
energies  of  the  bubble  and  any  heat  exchange  between  gas  and 
liquid.  We  further  assume  that  the  bubble  wall  remains  at  the 
liquid  temperature  T] ,  although  the  gas  and  vapor  mixture 
within  the  bubble  is  free  to  assume  any  temperature  dictated  by 
(6.23).  This  is  consistent  with  the  neglect  of  heat  transfer 
into  the  liquid.  The  partial  pressure  of  the  vapor,  Pv , 
within  the  bubble  is  taken  to  be  constant  and  equal  to  the 
equilibrium  vapor  pressure  at  Ti .  Thus  Pv  remains  constant 
due  to  the  assumption  that  the  bubble  wall  temperature  remains 
at  T^ .  Although  Pv  is  constant,  nv  will  change  as  the 
bubble  volume  changes  and  vapor  will  condense  or  liquid  vaporize 
at  the  bubble  wall.  This  results  in  a  form  of  heat  or  energy 
exchange  between  the  bubble  and  the  liquid  that  is  accounted  for 
by  the  term  nvhv  of  (6.23).  The  temperature  of  both  qas  and 
vapor  components  crossing  the  control  surface  in  either 
direction  is  taken  to  be  T^ ,  the  temperature  of  the  bubble 
wall . 


With  these  assumptions  and  the  previous  assumption  of  ideal 
gas  behavior  for  both  the  vapor  and  noncondensible  gas 
components,  we  can  write 


dU  =  c^ 


v,v 


d(n 
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v  ,g 
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Ui  t 

Here,  ci,j  i-s  the  mol ar  specific  heat  of  component  j  evaluated 
at  temperature  Tm.  The  notation  can  be  summarized  as: 

i  =  p,  specific  heat  at  constant  pressure, 

i  =  v,  specific  heat  at  constant  volume, 

j  =  g,  specific  heat  of  the  gas  component, 

j  =  v,  specific  heat  of  the  vapor  component, 

m  =  g,  evaluated  at  the  gas-vapor  mixture  temperature, 

Tg,  and 

m  =  1,  evaluated  at  the  bubble  wall  or  liquid  temperature, 

Ti. 

The  work  term  in  (6.23)  can  be  expressed  as  the  boundary  work  in 
moving  the  bubble  surface 


dw  =  (Pv  +  Pg)  dvb  . 


(6.26) 


Combining  (6.23-6.26)  and  rearranging,  we  obtain 


cg  Jr(nT)+cg  £_(nT) 
v,v  at  v  g  v,g  at  g  g 

=  -  (P  +  P  )  V.  +  (n  c1  +  n  c1  )  T.  . 

v  g  b  v  PfV  g  p  ^  g '  1 


(6.27) 


We  use  the  same  nondimensionalizations  as  before  and 


— m  m 

c .  .  =  c .  .  /R 

13  1,3  u 


(6.28) 
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We  then  write  (6.27)  in  dimensionless  form  and  obtain: 

cg  —  (n  T  )  +  cg  —  (n  T  ) 
v,v  dE  v  g'  v,g  g  9* 


\  <fv  +  V  +  <"v  c'p,v  +  "9  5p,g>  • 


(6.29) 


Recall  that  $  can  be  viewed  as  a  dimensionless  gas  constant  or 
as  the  ratio  of  PgQ  to  AP. 

Equation  (6.29)  together  with  the  ideal  gas  equations  of 
state  (6.18,  6.21)  are  the  equations  which  must  be  solved  for 
the  nonisothermal  case  in  place  of  (6.22)  for  the  isothermal 
case.  The  complete  set  of  simultaneous  equations  to  be  solved 
consists  of  the  bubble  dynamics  equations  (3.11),  the  boundary 
layer  approximation  for  gas  diffusion  subject  to  Henry's  Law 
(6.3)  and  conditions  (6.2,  6.3),  and  either  the  isothermal  ideal 
gas  change  of  state  equation  (6.22)  or  the  ideal  gas  equations 
of  state  (6.18,  6.21)  together  with  the  first  law  energy  balance 
(6.29).  This  set  of  equations  accounts  for  both  dynamic  and  gas 
diffusion  effects  on  bubble  behavior. 

6.3  Solution  of  the  Equations  Due  to  Mass  Diffusion 


We  now  describe  the  schemes  adopted  for  solution  of  the 
coupled  sets  of  equations  described  above  for  including  bubble 
dynamics  and  mass  diffusion  effects  in  bubble  behavior,  we 
first  approach  the  solution  of  the  energy  balance  (6.29)  by 
choosing  to  allow  it  to  determine  the  gas  partial  pressure  Pg 
and  its  rate  of  change  Pg .  With  this  scheme,  we  then  utilize 
(6.3,  6.13)  to  solve  for  ng  and  ng .  This  format  is  somewhat 
arbitrary.  However,  we  have  found  this  to  be  much  better 
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behaved  numerically  then  using  (6.3,  6.13)  to  determine  Pg 
differentiating  to  find  Pg  and  using  (6.29)  to  determine  n 
a  nd  ng . 

Combining  (6.21)  and  (6.29)  we  obtain 


cg  P  vw+cg  (P  V.+  V.P*  )  +  (P  +P  ) V,  —  6(n  c1  +n  c1  )=0  ,( 
v,v  v  b  v,gv  g  b  b  g'  v  g'  b  v  p,v  g  p,g  ,v 


where  the  condition  that  Pv  is  constant  is  employed.  We 
eliminate  nv  as  a  variable  by  relating  it  to  ng  using  (6.2 

"v  =  "g  5v/5g  ' 


"v  =  Pv  (n‘g/5g  -  "g  V^1  ‘ 


Thus  (6.30)  becomes 


4  [c9  P  V.  +  cg  (P  V,  +  V.  P  )  +  (P  +  P  )  V.  1 

6  L  v  j  V  V  D  y  '  g  q  q  g  '  '  v  ci  9  h  J 


g'  'b 


-1  v  •  .  -1 

c_  .  —  n„  +  c 
9 


P  P 

v _ g 


P'v  P 


PrV  g  p  2 


-1 

c  n 
Pf9  9 


=  0  .  ( 


The  relation  between  specific  heats  for  an  ideal  gas  in 
dimensional  form  is  given  by 

c  -  c  =  R  , 
p  v  u 

or  in  dimensionless  form 


c  =  c  - 


9 


6.30) 


1  ) 


6.31  ) 


6.32) 


1 


(6.33) 
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We  can  thus  express  all  values  of  cv  in  terms  of  the 
corresponding  Cp  quantities. 


Using  these  expressions,  equation  (6.32)  can  be  arranged  in  the 
form  of  a  differential  equation  for  Pg  : 


Pg  (A  +  B/P*)  +  Pg  D  -  E/Pg  +  F  =  0  ,  (6.34) 


where  A,...,F  are  defined  as: 


A  ■  '  ’>  VB 


*  (5?,g  -  "  - 


B  =  c  ^  n  P 
p,v  g  v 


D  *  cg  V  /B  *  cg  _  •  3  a  *  a  /  B  , 
p ,  g  b  p,g  o  o 


E  = 


-1 

c  n  P 
p,v  g  v 


F  = 


cg  P  V,  /  8  —  n 

p,v  v  b  p,g  9 


(6.35) 


Equation  (6,34)  is  solved  by  a  Runge-Kutta  procedure  when 
expressed  as 


E/P  -  D  P_  -  F 
9  9 

A  +  B/P* 


(6.36) 


where  A,...,F  are  taken  constant  and  known  for  a  given  time 
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-2  -  l  -g 

The  values  of  the  specific  heats,  Cpfg,  Cp^g  CpfV, 
and  CpfU,  are  taken  to  vary  only  with  temperature  (a 
"semi-perfect"  gas  assumption)  and  are  evaluated  using  the 
correlations  of  Van  Wylen  and  Sonntag  (1973,  p.  683). 

cp  =  (a  +  b  0P  +  d  0q  +  e  0r)/Ru  (6.37) 

where  6  =  T/100,  T  is  in  #K,  and  Cp  is  the  dimensionless 
specific  heat.  Values  of  the  constants  are: 


h2o 

a 

b 

d 

e 

P 

q 

n 

34.19 

-43.868 

19.778 

-0.88407 

0.25 

0.5 

i 

H 

8.9465 

0.0048044 

-42.679 

56.615 

1 . 5 

-1  .5 

-2 

n2 

9.3355 

-122.56 

256.38 

-196.08 

-1  .5 

-2 

-3 

These  correlations  are  stated  to  be  within  an  accuracy  of  less 
than  0.5  percent  for  temperatures  (in  degrees  K)  in  the  range 
300  <_  T  <_  3500.  For  temperatures  below  or  above  this  range,  the 
specific  heats  are  taken  as  being  constant  at  the  values  for 
T=300  or  T=3500  respectively. 

Air  is  treated  as  a  mixture  of  nitrogen  and  oxygen  such 
that  Cp(air)  =  0.21  Cp(02)  +  0.79cp(N2).  Since  the  liquid 
temperature  is  taken  to  be  fixed  for  a  given  problem,  Cp^g  and 
Cp,v  do  not  vary  and  are  evaluated  once  from  (6.37)  for  the 
specified  liquid  temperature.  The  specific  heats  dependent  on 
Tg  are  evaluated  each  time  (6.36)  is  solved  knowing  the 
dimensional  value  of  the  gas  temperature  Tg  =  Tg  •  T^ 
where  from  (6.21 ) 


Tracor  Hydronautics 


-95- 


T  =  T  •  T,  =  -2 —  •  T,  .  (6.38) 

S  ng 

The  specific  heats  are  evaluated  each  time  Equation  (6.36)  is 
solved  using  this  value  of  Tg  based  on  the  old  value  of  Pg . 

If  isothermal  conditions  are  assumed.  Equation  (6.22)  is  used  to 
solve  for  Pg. 

The  quantity  of  noncondensible  gas  in  the  bubble  at  any 
time  is  obtained  by  rewriting  equation  (6.13).  The  interval  of 
integration  in  (6.13)  can  be  split  into  two  intervals  to 
separate  out  the  current  time  step: 

1(0, t)  =  1(0, t-h)  +  I(t-h,t)  ,  (6.39) 

where  h  is  the  time  step  size  and 


h  na 

I  ( a  ,b)  =  r  - r- - -2 - ryj  dx  .  (6.40) 

a  [  f-x  ao  dy]  1/2 

The  second  integral  is  evaluated  to  avoid  the  singularity  in 
( t-x)-1/2 : 


( t-h,t)  =  j£_t 


n 


dx  = 


S02  (f-x',’/2 


-2  n 


5  ( t-x) 1/2 


x=t 


=  2  yK  n„/5^  .  (6.41) 
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where  the  tildes  denote  an  average  values  over  the  interval  from 
t-h  to  t.  We  approximate  these  average  values  as: 


ng  "  nq{t) 


(6.32) 


ao(t)  +  aQ(t-h) 


(6.43) 


Expression  (6.42)  for  rig  was  found  to  be  better  numerically 
than  the  mean  of  the  values  of  rig  at  t  and  t-h.  Then 


I ( t-h, t)  =  8 


[aQ(t)  +  aQ(t-h) ] 


(6.44) 


Equations  (6.13,  6.39,  6.44)  can  be  combined  to  provide  an 
expression  for  ng  (  t )  : 


ng  ( t )  = 


(aQ(t)  +  aQ ( t-h ) ] 


8  yrr 


i1  Cn{t)  ~  1(0, t-h)  ]  ,  (6.45) 


where 


u  t  D 

g 


1(0, t-h) 


/t-h 
>  o 


riU) 


dx 


(6.46) 
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Equation  (6.45)  is  used  to  calculate  rig(t).  The  value  of 
ng(t)  can  then  be  obtained  by  a  simple  trapezoidal  rule 
integration : 

n  (t)  =  n  (t-h)  +  h[n  (t)  +  n  (t-h)  ]/2 
y  y  y  y 

A  higher  order  approximation  could  also  be  easily  implemented. 
The  values  of  ng  and  ng  are  utilized  in  solving  (6.36)  or 
(6.22) . 

The  algorithm  employed  can  be  summarized  as  follows: 

1.  At  a  time  t* ,  all  variables  are  known. 

2.  At  time  t*+h,  the  bubble  dynamics  equations  are  solved 
for  Sq ,  aQ  with  all  other  parameters  evaluated  at  t=t*. 

Then,  the  differential  equations  for  all  other  components  of  the 
bubble  surface  equation  (see  Equation  3.10)  are  solved  with  the 
same  evaluation  of  the  parameters  at  the  preceding  time  step. 

3.  Equation  (6.45)  for  iig(t*+h)  is  solved  using  the 
values  of  aQ,  aQ  at  t*+h  and  the  value  of  Pg  at  t*.  It  is 
then  integrated  to  obtain  iig(t*+h). 

4.  Equation  (6.36)  or  (6.22)  is  solved  for  Pg  using  the 
values  of  aQ,  kQI  and  ng  at  t*+h  and  the  values  of  the 
other  parameters  at  t*.  Use  of  Henry's  law  yields  Cv(t*+h). 

5.  Iteration  is  performed  between  steps  3  and  4  above 
until  the  calculated  values  of  Pg  for  successive  iterations 

are  within  a  specified  difference  -  or  error  bound.  During  this 
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process,  the  "most  current"  values  of  the  variables  are  used  at 
each  step  of  the  iteration. 

6.  When  the  criterion  of  step  5  is  satisfied,  all 
variables  are  known  at  t=t*+h  and  the  calculation  returns  to 
step  1  with  all  variables  known  at  time  t=t*+h. 

This  process  is  continued  until  a  desired  time  is  reached. 

Figure  6.1  shows  a  flow  chart  of  the  program  GASCLD  based  on  the 
above  algorithm. 

6 . 4  Numerical  Examples  and  Interpretation 

A  series  of  numerical  cases  were  run  on  a  VAX  11/750 
computer  using  the  developed  program  GASCLD.  Only  a  few 
variables  were  investigated  in  order  to  demonstrate  the 
capabilities  of  the  code  and  to  gain  some  insight  into  the 
influence  of  gas  diffusion  on  bubble  cloud  dynamics.  Figures 
6-2  to  6-11  illustrate  the  results  obtained.  Bubble  growth  has 
been  considered  following  a  sudden  ambient  pressure  drop.  Two 
arbitrary  particular  values  of  the  initial  and  subsequent 
ambient  pressures  were  considered.  In  these  figures  comparisons 
are  made  between  the  behavior  of  an  isolated  dynamic-controlled 
bubble  (no  gas  diffusion),  an  isolated  bubble  with  gas 
diffusion,  and  a  multibubble  configuration  with  gas  diffusion. 

A  preliminary  investigation  of  the  influence  on  bubble  growth  of 
the  initial  concentration  of  dissolved  gas  in  the  liquid  and  of 
the  ratio  e  between  characteristic  bubble  size  and  characteris¬ 
tic  interbubble  distance  was  also  conducted. 

Figures  6-2  to  6-5  each  show  a  collection  of  cases  of 
normalized  bubble  radii,  r  =  r/rw  ,  versus  normalized  times, 


t 


.IUIUIL 
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t  =  t/r5  yp/(PQ  -  Pinf  )'•  1°  all  four  figures  the 

isolated  bubble  cases  are  presented  for  reference.  When  mass 
transfer  is  neglected  at  the  bubble-liquid  interface  a  poly¬ 
tropic  compression  law, 

P  Vk  =  (P  Vk)  at  t  =  0  ,  (6.48) 

g  g 


\ 

*m 

I 


has  been  adopted.  For  an  isolated  bubble  comparisons  are  made 
between  the  two  extreme  cases  of  no-gas  diffusion:  k  =  1 
(isothermal  behavior)  and  k  =1.4  (adiabatic  behavior),  and  the 
case  where  the  mass  transfer  and  the  energy  balance  equations 
are  solved.  For  both  cases  considered  here  the  adiabatic 
assumption  in  the  absence  of  gas  diffusion  gives  closer  results 
to  the  solution  of  the  diffusion  problem.  Based  on  the  perfect 
gas  law  (6.15),  this  result  is  to  be  expected  when  during  the 
bubble  growth  the  increase  in  the  number  of  moles  of  gas 
entering  the  bubble  overcomes  the  decrease  in  gas  temperature  in 
the  product  ( ng  Tg )  leading  therefore  PV  to  grow  with  V  (and 
not  remain  constant  as  k  =  1  presumes).  As  we  will  see  later  in 
Figure  6.10  the  temperature  drops  in  the  initial  phase  of  the 
growth  but  recovers  rapidly  to  maintain  its  value  later  in  the 
bubble  growth  history.  In  addition,  with  our  assumptions  of  no 
direct  heat  exchange  between  the  content  of  the  bubble  and  the 
liquid,  the  initial  phase  of  the  bubble  growth  is  adiabatic. 
Comparison  between  isothermal,  adiabatic  and  complete  solution 
results  show  that  deviation  between  the  three  cases  is  more 
significant  when  the  growth  is  pursued  for  a  much  longer 
nondimensional  time.  This  is  the  case  in  Figure  6.2  where,  the 
pressure  drop  being  less  violent  than  in  Figure  6.3,  the  growth 
rate  is  slower. 
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Comparisons  between  two  initial  gas  concentrations  can  be 
seen  in  Figure  6.2  As  expected,  for  an  isolated  bubble,  an 
increase  in  the  initial  gas  concentration,  CQ,  leads  to  an 
increase  in  the  bubble  radius  attained  at  any  time  mostly  due  to 
an  increase  in  the  amount  of  moles  transferred.  (As  we  will  see 
below  in  Figures  6.6  to  6.8,  this  change  in  CQ  affects  more 
significantly  the  effective  polytropic  coefficients.  In  a 
multibubble  system  interactions  become  very  significant  when  the 
bubble  sizes  approach  the  interbubble  distance.  This  is 
reflected  by  an  oscillation  in  the  bubble  radius  similar  to  that 
obtained  in  the  heat  transfer  problem.  Initially  during  the 
bubble  growth  the  modification  of  the  imposed  pressure  field  by 
the  presence  of  the  other  bubbles  in  the  cloud  reduces  the 
bubble  size  achieved.  At  any  instant  the  bubble  radius  is 
smaller  than  that  of  an  isolated  bubble.  Two  such  radii  are 
shown  in  the  figures  (R2  and  R  3)  and  correspond  to  the  sketch  in 
Figures  3.1  and  8.5.  R2  is  the  measurement  of  the  bubble  radius 
along  the  cloud  center  direction,  while  R 3  is  measured  in  the 
opposite  direction.  Figure  6.3  shows  the  variation  of  these  two 
dimensions  as  functions  of  time  for  a  set  of  values  of  the 
spacing  parameter  e.  It  is  obvious  from  this  figure  that  bubble 
interaction  and  deformation  increase  when  the  spacing  between 
bubbles  is  reduced.  The  presence  of  multiple  bubbles  introduces 
a  flattening  of  the  bubble  on  its  side  closer  to  the  cloud 
center.  This  side  has  a  slower  growth  than  the  other  side.  The 
computed  shrinking  and  oscillation  of  the  bubble  shape  in  the 
latter  stage  of  the  growth  needs  to  be  investigated  to  assure  it 
does  not  result  from  a  failure  of  the  asymptotic  approach.  This 
behavior,  seen  in  Figure  6.3  for  e  =0.2,  occurs  at  larger  times 
for  the  smaller  values  of  c.  All  the  computations  presented 
here  were  performed' with  a  logaLxthmic  increment  in  the  time 
steps.  Figures  6.4  and  6.5  present  some  of  the  results  of 
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Figures  6.2  and  6.3  with  logarithmic  time  scales.  A  smoother 
behavior  of  the  curves  can  be  seen  with  these  scales. 

Figures  6.6  to  6.10  show  for  various  conditions  the  time 
variations  of  the  effective  gas  polytropic  coefficients  when  gas 
diffusion  is  taken  into  account.  The  two  effective  coefficients 
used  in  these  figures  are  defined  as  follows.  The  total 
effective  polytropic  coefficient,  K  is  defined  base  on  the 
Relation  (6.48)  and  is  computed  at  each  time  step  using 

K  =  (Log  P  /P  ) / ( Log  V/V  )  .  (6.49) 

y  y  0 

The  local  effective  polytropic  coefficient,  K^,  is  based  on  a 
compression  law  which  applies  locally  between  two  computation 
steps 

(P  VKOi  =  (P  vK*)i  +  1  •  (6.50) 

It  is  therefore  computed  using  the  relation 

K  =  (Log  P  /P  ) / ( Log  V./V.  .)  .  (6.51) 

yi  yi  +  1 

Figures  6.6  to  6.9  show  that  the  total  effective  polytropic 
coefficient  drops  rapidly,  at  the  beginning  of  the  growth  from 
the  value  1.4  to  a  minimum  value  which  depends  on  the  initial 
gas  concentration,  CQ.  As  seen  in  the  figures  this  value  can 
become  significantly  smaller  than  one.  Since  at  these  later 
times  the  gas  temperature  becomes  practically  constant,  the 
lower  is  the  minimum  value  of  K  the  stronger  is  the  relationship 
between  the  number  of  moles  of  gas  diffusing  into  the  bubble  and 
the  bubble  volume  change.  This  can  be  reflected  by  the 
approximate  equation 
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(V/V  )1_K  =  n  /n  .  (6.52) 

y  yo 

derived  from  the  ideal  gas  law  with  Tg  assumed  constant.  This 
might  explain  why  the  minimum  value  of  K  is  smaller  when  CQ  is 
increased.  The  local  effective  coefficient  of  compression, 
decreases  much  more  significantly  than  K.  From  the  initial 
value  of  1.4,  decreases  practically  to  0  later  during  the 
growth  expressing  the  fact  that  Pg  and  Tg  become  practically 
constant  during  the  growth. 

The  variation  of  both  K  and  K4  are  much  more  sensitive  to 
the  initial  conditions,  (such  as  initial  gas  concentration, 

CQ,  and  initial  bubble  size,  Rq),  than  the  actual  bubble 
radius  history.  This  is  probably  due  to  the  fact  that  K  and 
Kjj  express  only  the  gas  thermodynamic  behavior  while  the 
bubble  radius  behavior  includes  gas  and  vapor  thermodynamic  and 
dynamic  phenoma  (pressure  drop).  These  results  are  well 
illustrated  in  Figure  6.9  where  the  influence  of  the  initial 
bubble  radius  is  seen.  No  difference  between  the  two  cases  of 
initial  bubble  radii  was  seen  in  the  radii  versus  time  curves. 
However,  the  effective  gas  polytropic  coefficients  are  very  much 
dependent  on  these  initial  conditions.  Figure  6.10  presents  the 
temperature  histories  for  these  two  initial  bubble  radii.  This 
figure  is  typical  of  all  temperature  versus  time  curves. 
Initially,  the  temperature  drops  significantly  and  attains  a 
minimum.  Later,  the  temperature  inside  the  bubble  recovers  to 
reach  the  initial  temperature  and  remains  practically  constant. 
The  minimum  temperature  achieved  and  the  time  to  achieve  this 
minimum  are  functions  of  the  initial  conditions.  This  result  is 
related  to  the  variations  of  K  and  Kt.  The  influence  of 
temperature  change  on  the  effective  coefficients  K  can  be  seen 
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in  Figure  6.11  where  the  solution  based  on  the  energy  equation 
is  compared  with  that  where  the  temperature  was  imposed  to 
remain  constant. 

The  above  results  should  be  extended  to  actual  pressure 
fields  where  the  pressure  variations  are  more  realistic  than  a 
sudden  pressure  drop  and,  more  importantly,  to  longer 
computation  times  where  gas  diffusion  effects  become  more 
significant.  For  such  more  practical  situations  the  influence 
of  the  magnitude  of  the  pressure  drop,  P0  -  P^nf,  should  be 
investigated  carefully. 
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7.0  CONTINUUM  MEDIUM  APPROACH 


7.1  Introduction 


One  major  assumption  of  the  approach  presented  above  is 
that  in  first  approximation  the  imposed  ambient  pressure  is 
assumed  to  be  instantaneously  transmitted  to  the  vicinity  of 
each  bubble  in  the  cloud.  Therefore,  both  the  compressibility 
of  the  bubble  medium  and  the  influence  of  the  liquid 
motion  generated  by  the  other  bubbles  on  the  dynamics  of  the 
bubble  considered  were  neglected  in  the  first  order 
approximation.  This  limits  the  validity  of  the  study  to  very 
low  void  fractions.  The  incompressibility  assumption  is  valid 
as  long  as  the  fluid  velocity  does  not  approach  the  speed  of 
sound.  For  single  bubble  dynamics,  this  does  not  usually  happen 
until  the  final  phase  of  the  collapse.  Here,  however,  two 
factors  contribute  to  limit  the  validity  of  the  assumption. 
First,  the  rate  of  implosion  is  higher  for  a  cloud  than  for  a 
single  bubble,  and  second,  more  importantly,  the  velocity  of 
sound  drops  considerably  when  the  void  fraction  increases.  This 
underlines  the  need  to  account  for  the  behavior  of  the  cloud  as 
a  whole  in  order  to  determine  a  more  accurate  value  of  the  local 
pressure  driving  the  collapse  of  the  individual  bubbles.  In 
addition,  this  would  have  the  advantage  of  limiting,  for  the 
following  orders  of  approximations,  the  number  of  bubbles 
directly  influencing  the  considered  one  through  a  time  delay  of 
the  propagation  of  the  information  from  one  bubble  to  another. 
Indeed,  the  asmyptotic  theory  shows  that  the  effective  parameter 
of  the  expansions  is  ec  i ,  (where  clr  defined  by  (3.13),  is  a 
direct  function  of  the  number  of  bubbles),  rather  than 
c  =  r b0/*0*  Introducing  a  motion  equation  for  the  bubbly 
medium  would  limit  the  number  of  influencing  bubbles  to  those  in 
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the  direct  vicinity  of  the  considered  one,  through  a  time  delay 
of  the  propagation  of  the  information  from  one  bubble  to 
another.  In  summary,  if  we  account  for  a  motion  equation  in  the 
cloud  medium,  the  first  order  approximation  of  the  preceding 
approach  becomes  more  accurate,  and  as  a  consequence  the 
following  corrections  will  be  smaller,  making  the  approach  valid 
for  higher  void  fractions,  a. 

7 . 2  Volume  Average  Method 

Basically  the  classic  methods  used  to  describe  a  two-phase 
medium  are  not  much  different  from  the  singular  perturbation 
method  presented  above.  The  final  description  deals  only  with 
the  macroscale  of  the  cloud. 

The  description  of  the  macroscale  of  the  cloud  can  be 
obtained  by  averaging  the  various  physical  quantities  defined  in 
the  microscale.  The  two-phase  medium  is  assumed  to  be 
constituted  of  "particles"  containing  the  host  liquid  and  few 
bubbles.  This  "particle"  is  small  enough  to  be  able  to 
distinguish  the  gaseous  and  liquid  constituents,  but  large 
enough  to  enable  one  to  define  significant  volume  average 
quantities  in  the  two-phase  continuum.  Therefore,  each 
"particle"  appears  in  the  macroscale  as  a  fluid  point  M  allotted 
various  physical  and  kinematic  properties:  a(M,t)  is  the  local 
void  fraction,  pm(M,t)  is  the  local  medium  density  ,  Um(M,t) 
is  the  velocity,  and  Pm(M,t)  the  pressure,  ...,etc.  In  such  a 
volume  averaging  description,  if  Vp  is  the  volume  of  the 
particle,  X(M,t)  the  considered  average  quantity  and  x(m,t)  its 
local  value  in  the  microscale,  we  have  the  following  definition: 
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X  ( M ,  t )  =  i  /  x (m , t )  d  V  (7.1  ) 

P  V 

P 

The  density  of  the  medium  is  therefore  defined  by  the 
relation: 

Pm(M,t)  =  [l-a(M,t)]  +  pg  a(M,t)  .  (7.2) 

The  liquid  is  assumed  to  be  incompressible  and  p4  constant. 

The  void  fraction,  a(M,t),  is  defined  as  the  relative  volume  of 
gas  in  the  particle.  Usually  pg a  is  neglected,  and  the 
density  of  the  medium  is  written: 

Pm(M,t)  =  Pjl  [1  -  a(M , t )  ]  .  (7.3) 


If  U£(M,t)  is  the  average  velocity  of  the  liquid  in  the 
particle  and  Ug(M,t)  the  average  velocity  of  the  gas,  we 
obtain  comparable  results  to  (7.2)  and  (7.3): 


pm  “m  =  p*“  £  <1- 

a)  +  Pg  Ug  a  , 

(7.4) 

pm  “m  *  p£z  - 

a)  » 

(7.5) 

and  combining  with  (7.3), 

"  * 

(7.6) 

The  continuity  equation  is  obtained  by  writing  the  mass 
conservation  of  a  volume  of  the  bubble  medium  followed  during 
its  motion.  Using  the  average  quantities  defined  above  we  can 
write: 
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Here  the  material  derivative  pertains  to  the  medium  velocity 
Um  or,  with  our  assumptions,  to  U*  because  of  Equation 
(7.6): 

d/dt  =  3/  3t  +  U£  V  .  (7.8) 


Since  Equation  (7.7)  is  valid  for  any  volume  V,  we  obtain  the 
general  equation: 

TT  +  7  •  <  W  '  0  '  l7'9» 

A  similar  equation  can  be  written  concerning  the  number  of 
bubbles,  n(M,t).  Neglecting  any  complete  bubble  disappearance 
or  sudden  generation,  as  well  as  bubble  splitting  and 
coalescence,  we  can  write: 

^  +  n  7  •  =  0  ,  (7.10) 

the  material  derivative,  D/Dt,  being  defined  as: 

D/Dt  =  3/ 3t  +  Ug  V#  .  (7.11) 

The  momentum  equation  of  the  bubbly  medium  can  be  obtained 
in  the  same  manner  by  using  the  momentum  equations  of  both 
constituents  in  the  microscale  and  integrating  over  the 
"particle"  volume  Vp.  If  we  neglect  the  viscous  forces,  this 
can  be  written: 

d .  U. 

J  [  pi  “at1  +  ^  ]  d  V  =  0  ,  (7.12) 
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the  index  i  designating  the  liquid  or  the  gaseous  phase 
depending  on  the  position  of  the  element  of  volume  dV  in  the 
microscale.  If  we  account  for  the  incompressibility  of  the 
liquid  this  equation  becomes: 


°t  "  •  a) 


dt 


+  p  a 

g 


Dt 


-  Pg  /  Ug  dV  /  =  0 

VP  VP  (7.13) 


If  we  neglect  the  gas  contribution  to  the  momentum,  and  we 
account  for  (7.6)  we  obtain  the  following  approximate  classical 
momentum  equation: 


(1  -  ex) 


dU 
— m 

dt 


+  V  P  =  0 
m 


(7.14) 


where  it  is  assumed  that 

/  Ap.  d  V  =  /  p.  n  ds  »  Vp  •  V  Pm  .  (7.15) 

VP  AP 

7 . 3  Bubble-Liquid  Relative  Motion 

The  only  equation  left  is  that  giving  the  bubble 
translation  velocity,  Ug ,  which  reflects  the  interaction 
between  the  two  phases  of  the  bubbly  medium.  The  study  of  this 
equation  is  a  whole  subject  of  research  in  itself.  Several 
contributions  exist  which  have  dealt  with  more  and  more 
complicated  situations.  When  viscous  drag  is  neglected,  a  very 
interesting  general  expression  for  the  motion  of  a  deformable 
bubble  in  a  nonuniform  potential  flow  was  derived  by  Landweber 
and  Miloh  (1980).  If  we  admit,  however,  that  the  liquid  flow 
around  an  isolated  bubble  is  linearly  accelerated,  and  that  the 
bubble  remains  in  first  approximation  spherical,  we  can  write, 
neglecting  the  bubble  mass,  a  simpler  equation  as  follows: 
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DU„  DU , 

~2.  -  3  - - 

Dt  Dt 


V 


(7.16) 


In  this  equation,  the  virtual  mass  of  the  bubble  is  considered 
to  be  2/3na3Qp£  and  the  material  derivative  is  related  to 
the  bubble  velocity  as  discussed  by  Prosperetti  and  Van 
Wijngaarden  (1976). 

When  other  bubbles  are  present  in  the  flow,  corrections  are 
to  be  introduced  in  this  expression,  following  Landweber's 
calculations.  Van  Wijngaarden  (1976)  performed  similar 
corrections  for  a  rigid  sphere  and  obtained  the  expression: 


d 

“Ht 


1 

2 


(  1  +  2.78  at)  (U^ 


‘ 


]»(!-«) 


~St  V-l 


(7.17) 


A  relationship  similar  to  (7.16)  or  (7.17)  was  implemented 
in  the  program  CLDMAIN2  to  account  for  gravity  effects. 


7.4  Micromorphic  Continuum  Description 

In  classical  continuum  mechanics  the  fluid  is  described 
geometrically  by  a  field  point  M  and  kinematically  by  a  velocity 
field  U ( M ) .  The  averaging  approach  of  the  cloud  medium,  as 
described  in  the  preceding  paragraph,  is  in  this  sense 
classical.  However,  when  a  medium  contains  microstructure,  as 
is  the  case  for  a  bubbly  medium,  a  more  refined  description  can 
be  obtained  by  assigning  to  M,  in  addition  to  the  macroscale 
velocity,  U(M) ,  other  quantities  which  reflect  the  microscale 
behavior  in  the  "particle".  In  a  first  gradient  theory,  in 
addition  to  the  velocity  field,  U(M) ,  a  field  of  the  gradients 
of  relative  velocities  in  the  microscale  scale,  X,  is  added 
which  defines  kinematically  the  medium.  The  description  can  be 
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further  refined  by  using  higher  order  gradient  theories. 

Germain  (1973)  considered  such  approaches  and,  using  the  method 
of  virtual  power,  was  able  to  derive  the  equation  of  motion  of 
the  continuum  medium  accounting  for  the  macrostresses,  _cr,  and 
the  microstresses,  £. 

In  a  first  gradient  theory  the  velocity  in  the  microscale 
can  be  written  as 

U'(m)  =  U  ( M)  +  ^M)  •  Mm  ,  (7.18) 

where  m  is  a  point  in  the  microscale  (see  Figure  7.1). 
Consequently  the  acceleration,  _T*  ,  of  m  is  derived  and,  by 
equating  at  dynamical  equilibrium  the  virtual  power  of  all  the 
internal  and  external  forces  acting  on  the  considered  particle 
(volume  Vp)  to  the  material  derivative  of  the  virtual  power  of 
mass  velocity  of  Vp,  one  obtains  a  dynamical  equation  of  the 
medium  relating  S,  a,  and  JT *  . 

To  define  X  we  consider  the  motion  on  a  scale  which  is  of 
the  same  order  as  the  microstructure.  To  do  so  for  a  bubble 
cloud,  let  us  divide  the  cloud  medium  into  fluid  "cells"  each 
enclosing  an  isolated  bubble.  In  addition,  we  assume  for 
simplicity  that  the  bubble  center  of  mass  and  the  "cell"  center 
of  mass  coincide  at  the  considered  time.  Let  U(M)  be  the 
velocity  in  M  induced  by  the  rest  of  the  cloud  in  absence  of  the 
bubble,  and  V(B)  the  velocity  of  the  bubble  center,  B.  U(M) 
would  be  the  value  of  the  velocity  field  assigned  to  M  in  a 
classical  fluid  mechanics  description. 

The  bubble  radius  is  aQ  and  its  variations  with  time  are 
given  by  the  Rayleigh  Plesset  equation.  This  radial  motion  of 
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the  bubble  surface  induces  at  a  point  m  of  the  cell  (Figure  7.1) 
a  velocity  of  value  (ur  •  er'),  where  er’  is  the  unit 
vector  of  the  direction  Mm.  The  total  velocity  u' ,  at  m  is: 


u' (m)  =  U ( M)  + 


•  2 

Vo 

r'2 


e  ' 
— r 


v  [- 


2r 


2  - 


(U(M)  -  V ( B ) ) 


±r’  ]  + 

(7.19) 


where  r'  is  the  distance  between  M  and  m.  The  second  term  in 
this  expression  is  a  source  term  due  to  the  spherical  bubble 
oscillations,  while  the  last  term  is  a  dipole  due  to  the  slip 
velocity  between  the  spherical  bubble  and  the  fluid,  and  could 
include  first  order  corrections  of  the  bubble  shape.  For 
further  corrections  for  nonsphericity  of  the  bubble,  other  terms 
(singularities  of  higher  orders)  have  to  be  included.  By 
differentiating  (7.19)  with  respect  to  time  and  space  one  can 
define  an  acceleration  vector,  r' ,  and  a  strain  rate  tensor, 

D'.  Following  Germain's  approach,  and  using  the  principle  of 
virtual  powers,  one  could  then  derive  an  equation  of  motion  of 
the  cloud  medium.  We  decided  instead  to  start  with  a  first 
gradient  theory  and  replace  (7.20)  by  its  Taylor  expansion.  We 
follow  in  doing  so  the  first  calculations  done  by  Michelet 
(1980)  in  his  graduate  thesis. 


The  basic  approximation  used  in  this  linearization  approach 
is  based  on  the  fact  that  Equation  (7.19)  is  only  valid  in  the 
liquid  portion  of  the  "cell"  (r  >  ao).  It  seems  therefore 
logical  to  write  the  velocity  at  m,  close  to  the  bubble 
boundary,  as  a  Taylor  expansion  of  the  value  of  u'  computed  on  a 
point  of  the  bubble  surface,  S,  (figure  7.1).  This  has  the 
advantage  of  eliminating  the  singularity  of  (7.19)  for  r'  =  0. 
The  obtained  expression  for  u'(m)  is  then: 
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u'(m)  =  U  { M)  +  [3aQ  +  4Vfc  cos  0  -  r'  ^2  ~  +  3  ~  cosej]  er 


+  [2  Vt  sin8  -  r'  Vfc  sin0]  e( 


(7.20) 


where  Vt  =  V(B)  -  U(M) | . 

When  Vt  is  not  accounted  for,  the  expression  of  u' (m)  reduces 
to  a  form  comparable  to  (7.18),  which  is  much  easier  to 
interpret  than  equation  (7.20).  In  that  case  we  obtain: 

u  '  (M‘  )  =  U(M)  +  ^  •  Mm  +  sl  *  e  (7.21  ) 


where  X  and  oi  are  both  tensors  assigned  to  M  and  defined  as: 

a 

-2  ~  L  f  £=  =  3  aD  L  '  (7.22) 

I.  is  the  spherical  unit  tensor.  We  notice  that  in  comparison  to 
(7.18),  which  describes  a  first  gradient  homogeneous 
deformation,  in  (7.21)  there  is  in  addition  to  the  gradient 
tensor,  X,  a  tensor  a  reflecting  the  presence  of  a  source  in  the 
cell.  Equation  (7.20)  reflects  in  addition  to  this  the  presence 
of  a  dipole.  It  could  be  written  as 

u’(m')  =  U ( M )  +  jL  •  Mm  +  X_'  •  |  Mm  |  !z  +  j.  •  er  +  •  ez 

+  (6Vfc  cose  -  r'  Vt  cos  9)  er  ,  (7.23) 


where  ez  is  the  unit  vector  of  the  direction  of  U  and  V;  X' 
and  a'  play  the  same  role  as  X  and  a  but  are  applied  just  to  the 
direction  of  the  translation.  The  last  inhomogeneous  ter.i  is 
more  difficult  to  put  in  simple  form. 
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From  the  expression  (7.20)  we  can  now  compute  the 
acceleration,  then  apply  the  principle  of  virtual  power,  to 
obtain  the  equation  of  motion.  Here  again  in  the  absence  of 
translation  velocity  ,  the  results  are  simpler  to 
interpret.  In  the  absence  of  viscous  effects  these  results  can 
be  written  as  follows: 


m 


dU 

dt 


+  3  K  (h 


(7.24) 


where  K  depends  unfortunately  on  the  cell  geometry 


p  _  K  =  /  p-  e  d  V 

m  -  Jv  l  -r 


(7.25) 


c 

If  the  cell  and  the  bubble  are  symmetrical  with  regard  to  the 
center  of  mass  M,  then  K  =  0,  and  (7.25)  reduces  to  the 
classical  equation  of  motion  (7.14).  Although  it  is  unfortunate 
that  the  cell  shape  seems  to  play  a  role  in  the  model,  K  might 
rather  reflect  an  effect  of  the  nonsphericity  of  the  bubble. 


When  Vt  is  taken  into  account  a  whole  series  of  "inertia" 
integrals  like  (7.25)  appear  in  the  calculations.  In  order  to 
see  what  such  a  model  might  indicate  we  considered  the  case  of  a 
spherical  bubble  in  a  spherical  cell.  In  this  case  the  motion 
equation  becomes: 


0(  a)  ] ( V  -  U)  }  =  -7p. 


(7.26) 


Here,  R  is  the  radius  of  the  cell,  and  if  we  write 
R  *  a0a-^/3,  we  have  the  result: 
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dU 


m  ldt 


+  3 


r  -1/3  ,1  , 

[a  (t  +  a 


’/3  -  2a2/3  + 


] (V  -  U)  }  =  -  Vp 
(7.27) 


This  surprising  result  (dependence  on  a”^/3)  might  be  compared 
with  that  obtained  for  the  apparent  viscosity  of  a  bubbly  flow, 
which  is  4 u/3  •  a-1  (Batchelor  (1967),  Van  Wijngaarden  (1972)). 
In  addition,  due  to  the  linearization  of  the  velocity  field 
(first  gradient  theory)  this  model  loses  it  validity  for  low 
a '  s . 


7 • 5  Example  of  a  Spherically  Symmetrical  Cloud 

Let  us  consider  a  finite-sized  spherical  cloud  of  bubbles 
and  define  its  radius,  R(t),  at  time  t,  as  the  position  of  the 
last  outer  shell  of  bubbles.  The  space  is  therefore  divided 
into  two  regions.  For  r  >  R(t),  the  medium  is  an  incompressible 
liquid  of  density  p£,  while  the  interior  of  the  sphere, 
r<R(t),  is  filled  with  a  two-phase  medium  which  can  be  defined 
as  in  Section  7.3.  We  define  at  a  point  M(r),  a  radial  liquid 
velocity  Ujj(r,t)  and  a  radial  bubble  translation  velocity 
Ug(r,t).  Similarly,  we  define  a  local  void  fraction  a(r,t), 
density  pm(r,t),  bubble  radius  a0(r,t),  number  density  n(r,t), 
and  medium  velocity  um(r,t).  The  matching  between  the  two 
media  states  that  at  r  =  R(t)  there  is  a  continuity  of 
velocities  and  pressures: 


R( t )  =  u  ( r  , t )  ,  (7.28) 

P  .(R,t)  =  lim  p’ (R,  r' ,  t)  ,  (7.29) 
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where  r'  is  the  distance  in  the  microscale  between  a  babble 
center  and  a  cell  field  point.  The  continuity  and  momentum 
equations  in  the  liquid  medium  (r>R(t))  are  easy  to  solve  and 
give,  after  neglecting  viscous  effects: 


u£(  r  )  =  Vg/4  nr  2 


(7.30) 


(7.31  ) 


Vg  is  the  total  volume  of  the  bubbles  in  the  cloud, 

V  =  4i  f  a  r2  dr  .  (7.32) 

y  o 


Inside  the  bubbly  medium,  due  to  the  spherical  symmetry, 
the  continuity  equation  also  gives 


with 


u£(r,t) 


V 


(r  ) 


4  nr 


r  <  R(t) 


V(r)  =  4  n  /r  a{  x  ,  t )  x2  dx 
o 


(7.33) 

(7.34) 


If  we  are  interested  in  the  problem  of  the  collapse  of  the  cloud 
under  an  imposed  ambient  pressure  variation,  P^t),  (20)  can 
be  integrated  between  the  cloud  radius  and  infinity  to  give: 


-PJt)  +  P£(R) 


(7.35) 


Using  (7.29),  P£(R)  can  be  related  to  the  behavior  of  any 
individual  bubble  of  radius  a0  in  the  last  outer  shell  of  the 
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cloud,  using  the  Rayleigh-Plesset  equation.  Equation  (7.35) 
becomes : 


(7.36) 


The  cloud  radius  motion  can  be  obtained  by  using  an  equation  of 
the  bubble  motion,  for  example,  (7.16).  This  gives  the 
following  second  relation  between  R,  a0,  and  Vg: 


•  • 

R  + 


V9  2 
2R  3 


(7.37) 


A  third  equation,  in  addition  to  (7.36)  and  (7.37),  is 

needed  to  solve  for  R,  a0,  and  Vg .  Without  an  assumption  on  a 

proportionality  between  Vg(t)  and  a0(t)  or  without  penetrating 

the  cloud  and  solving  for  all  a0(r,t)  to  determine  Vg ,  there 

is  no  hope  of  solving  the  problem.  We  do  not  think  the 

proportionality  assumption  is  generally  justifiable  even  if  at 

t  =  0  all  bubbles  in  the  cloud  have  the  same  size,  since  P(r,t) 

would  not  generally  be  the  same  for  any  location  r  at  a 

subsequent  time.  However,  an  assumption  similar  to  that  done  by 

•  ■)  • 

Morch  (1980)  could  be  done  stating  the  Vg  =  n  na  a  where  n  is 
the  number  of  bubbles  in  the  outer  shell.  The  need  in  general 
to  solve  the  whole  problem  is  to  be  expected  and  is  very 
important  because  it  shows  that  defining  the  cloud  by  just  one 
parameter,  as  a  unique  void  fraction,  is  not  sufficient  to 
describe  its  dynamics.  Number  and  bubble  size  distribution  are 
other  important  variables  to  consider. 
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The  problem  presented  above  can  be  solved  by  combining  the 
multibubble  approach  with  the  global  descriptions  presented  in 
Sections  3.1  and  3.3.  The  approach  of  d'Agostino  and  Brennen 
(1983)  is  very  close  to  that  presented  in  this  section  with  the 
additional  assumption  of  small  radius  and  pressure  oscillations 
and  the  neglect  of  bubble-liquid  relative  motion  and  bubble 
interactions . 
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8.0  CLOUD  CAVITATION  THROUGH  A  SIMULATED  MULTIBUBBLE  SYSTEM 
8 . 1  Description  of  Experiment 

Since  it  is  difficult  to  measure  the  quantities  involved 
within  an  actual  cavitation  cloud  produced,  for  example,  .n  a 
cavitation  tunnel,  we  have  resorted  to  the  same  techniques  used 
to  study  single  bubble  collapse,  namely  spark  generated 
bubbles.  The  fast  discharge  of  high  voltage  condensors  across 
two  submerged  electrodes  produces  a  bubble  which,  once  attaining 
its  maximum  volume,  is  not  in  equilibrium  with  the  ambient 
pressure.  The  growth,  rapid  collapse,  and  motion  of  the  bubble 
can  be  recorded  by  high-speed  photography.  In  addition,  the 
pressure  fluctuations  caused  by  the  bubble's  dynamics  can  be 
monitored  by  means  of  a  pressure  transducer.  The  idea  we  have 
used  based  on  this  principle  consists  of  placing  a  system  of 
nearby  solid  walls  which,  through  reflections,  generates  a 
fictitious  cloud  of  bubbles  composed  of  the  actual  bubble  and 
its  images.  We  know  from  the  method  of  images  in  potential 
flows,  that  solid  walls  behave  schematically  as  mirrors.  For 
this  study,  wall  structures  which  depict  symmetrical 
configurations  of  N  =  4,  8,  12,  and  N  >  12  bubbles  were  used. 

A  vacuum/pressure  tank  was  built  containing  the  electrodes, 
transducer,  hydrophone,  and  one  of  the  various  multiple  wall 
structures.  The  location  of  the  spark  gap  was  varied  in  the 
test  tank  in  order  to  control  6,  the  ratio  of  the  maximum 
characteristic  radius  of  the  bubble  produced,  Rq,  to  the 
distance  between  the  electrode  tips  and  the  apex  of  the  wall 
structure,  ds.  6  is  directly  related  to  e,  the  ratio  between 
bubble  characteristic  radius  and  the  bubble  spacing.  This 
relation  is  a  simple  factor  of  proportionality  which  is  geometry 
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dependent.  The  spark-generating  system  was  capable  of  obtaining 
a  voltage  charge  up  to  16  kV  on  a  0.95  pfd  capacitor.  For  the 
purpose  of  these  tests,  a  charge  voltage  of  12  kV  was 
consistently  used.  The  system  was  made  more  versatile  through 
the  addition  of  an  electronic  circuit  which  allows  it  to  fire  at 
a  controlled  repetition  rate  which  could  be  of  a  high 
frequency.  The  relative  locations  between  the  components  were 
varied  independently  so  as  to  study  the  effect  of  changing  the 
distance  between  the  bubble  and  the  walls. 

The  signals  from  the  transducers  and/or  the  hydrophone  were 
captured  on  a  Nicolet  2090  digital  oscilloscope  and  then  stored 
on  floppy  disks  for  later  analysis.  This  oscilloscope  has  a 
maximum  digitizing  rate  of  50  ns  per  point  and  a  bandwidth  of  up 
to  7  MHz  with  a  vertical  resolution  of  0.4  percent.  The 
oscilloscope  was  triggered  either  from  the  initial  pressure 
pulse  caused  by  the  spark  or  from  a  signal  from  a  variable  time 
delay  circuit  if  the  collapse  of  the  bubble  was  to  be 
investigated  using  an  expanded  time  scale.  A  Digital  Equipment 
Corp.  VAX  11/750  computer  was  interfaced  with  the  oscilloscope 
allowing  the  digitized  pressure- time  data  to  be  analyzed  with  a 
Fast  Fourier  Transform.  Figure  8.1  is  a  photograph  of  the 
overall  spark  bubble  test  facility  showing  the  spark  generator, 
test  tank,  and  data  acquisition  equipment.  Figure  8.2  is  a 
photograph  showing  the  test  tank,  the  electrodes,  transducer, 
and  the  multiple  walls  structure,  in  this  case  a  cone.  Figure 
8.3  shows  a  representative  curve  of  the  pressure  time  signal  for 
a  spark  induced  bubble  and  Figure  8.4  the  Fourier  transform  of 
the  same  data.  For  this  particular  case,  the  electrodes  were 
located  2.25  in.  from  the  top  of  the  cone.  The  ambient  pressure 
inside  the  tank  was  13.1  kPa. 
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In  this  study,  we  tested  a  two-sided  90  degree  corner  (two 
solid  plates  perpendicular  to  each  other),  a  three-sided  90 
degree  corner  (three  solid  plates  perpendicular  to  each  other), 
a  pentagonal  base  pyramid  (five  equi-sized  triangles  held 
together),  and  a  cone  with  a  60  degree  angle  at  its  apex.  A 
series  of  tests  was  run  varying  the  distance  from  the  apex  of 
each  of  the  wall  structures  to  the  electrodes.  For  all  tests, 
the  transducer  was  fixed  at  6  inches  below  the  electrodes,  and 
the  ambient  pressure  was  approximately  14  kPa.  High-speed 
movies  were  taken  of  the  bubble  growth  and  collapse  for  each 
test,  and  the  bubble  radius  history  was  measured  from  the  films 

8.2  Discussion  of  Results 


To  analyze  the  high-speed  films  of  the  spark  generated 
bubbles,  the  model  shown  in  Figure  8.5  was  used  consistently, 
and  shall  also  be  used  for  the  results  presented  in  this 
report.  It  shows  the  selected  characteristic  dimensions  of  the 
bubble  shape  measured  relative  to  the  initial  bubble  center,  B, 
and  the  apex  of  the  wall  structure,  point  O.  Point  O  also 
represents  the  center  of  the  fictitious  symmetrical  cloud 
configurations.  Variation  of  R2  with  time  represents  a 
measurement  of  the  re-entering  jet  advancement.  From  the 
high-speed  films,  it  was  observed  that  for  all  tests, 
theoretically  predicted,  this  re-entering  jet  formed  and  was 
always  oriented  in  the  direction  of  the  apex,  or  cloud  center, 
point  0. 


The  details  of  this  jet  formation  (time  of  formation  and 
speed  of  advancement)  differs  significantly  between  the  various 
test  cases.  Figures  8. 6-8. 9  show  the  time  histories  of  R2 
measured  for  several  values  of  0  in  configurations  of  N  =  4,  8, 
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12,  and  N  >  12  bubbles.  These  configurations  correspond  to  the 
wall  structures  of  the  2-D  corner,  the  3-D  corner,  the 
pentagonal  pyramid,  and  the  60  degree  cone,  respectively.  All 
the  curves  have  been  normalized  using  the  following  procedure 
which  tends  to  minimize  the  influence  of  spark  energy  variations 
and  ambient  pressure  variations  between  one  case  and  another. 

All  bubble  dimensions  were  normalized  with  an  average  maximum 
bubble  radius  Rq, 


Ro  =  (2Ri  +  R2  +  R  3 )/4 . 

Time  was  also  normalized  using  Rq,  the  measured  absolute 
pressure  in  the  tank,  Pa,  and  the  vapor  pressure,  Pv .  The 
characteristic  time  used  for  this  purpose  is  the  Rayleigh  time, 

tR' 


where  p  is  the  liquid  density. 

Close  examination  of  these  plots  suggests  that  increasing 
the  number  of  bubbles  increases  the  strength  of  the  re-entering 
jet,  ( R 2  becomes  more  negative).  While  the  results  from  the 
simulated  4  bubble  configuration  show  that  the  jet  did  not 
penetrate  past  the  original  center  of  the  bubble,  B,  at  all,  the 
results  from  the  simulated  12  bubble  configuration  and  the 
circular  cone  tests  show  a  penetration  past  this  point  of  more 
than  40%  of  the  maximum  radius  in  some  cases.  Since  the 
strength  of  the  re-entering  jet  gives  an  indication  of  the 
degree  of  erosion  from  cavitation,  these  results  verify  the 
theoretical  predictions  about  the  increased  cavitation  erosion, 
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( increased  strength  of  the  collapse),  from  a  cloud  as  compared 
to  a  single  bubble. 

Another  trend  which  is  noticeable  from  these  results  is 
that  the  jet  penetrates  further  inward  for  higher  values  of  B 
(or  e) .  This  trend  is  shown  in  Figures  8. 6-8. 9  for  all  of  the 
bubble  configurations,  but  becomes  much  more  noticeable  for  the 
higher  density  configurations  (pyramid  and  cone  tests)  and  the 
closer  interacting  bubbles.  Since  higher  0's  correspond 
theoretically  to  higher  void  fractions,  this  result  agrees  with 
the  theory  in  that  stronger  collapses  are  occurring  for  these 
cases.  This  is  also  shown,  perhaps  more  clearly,  in  Figure 
8.10,  which  presents  the  variation  of  the  minimum  value  of  R2 
with  6.  One  can  also  see  that  for  the  N=4  and  N=8  bubble  cases, 
a  strong  minimum  also  occurs  at  0.09  <  0  <  0.11.  This  minimum 
is  noticed  slightly  for  the  N=12  bubble  case  and  not  at  all  for 
N  >  12  bubbles.  This  result  which  proved  to  be  hard  to  check 
with  certainty  with  our  testing  facility  indicates  that  there 
exist  optimum  values  of  0  (or  e)  for  maximum  collective  effect 
(increased  collapse  intensity,  erosion,  or  noise  capability). 

Figures  8.11-8.14  show  the  radius  histories  for  R3 ,  the 
measurement  of  the  bubble  radius  along  a  direction  toward  the 
center  of  the  cloud.  From  these  results,  one  can  see  that  this 
measurement  is  closely  tied  to  R2.  After  the  first  collapse,  R3 
is  seen  to  grow  with  increased  0,  and  the  bubble  moves  closer  to 
the  center  of  the  cloud.  Again,  the  net  result  is  that  the  void 
fraction  for  the  simulated  cloud  increases,  subsequently 
increasing  the  strength  of  the  collapse  as  long  as  the  generated 
bubble  does  not  touch  the  walls. 


I 
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Figure  8.15  shows  the  period  of  bubble  oscillation  with 
respect  to  8  for  the  various  symmetrical  configurations.  It 
suggests  a  lengthening  effect  in  the  bubble  period  for  small  and 
large  values  of  8,  with  a  minimum  period  occuring  in  the  region 
.07  <  8  <  .13  for  all  the  cases  considered.  One  can  see  that 
the  strongest  minimum  occured  during  the  tests  run  with  the  60° 
circular  cone.  The  data  from  this  figure  and  that  of  Figure 
8.10  suggest  an  optimum  value  for  8  for  maximum  noise  and 
erosion  within  the  range  considered  for  the  test.  This  optimum 
seems  to  be  at  or  near  8  =  0.1  for  all  of  the  bubble 
configurations.  This  result  is  consistent  with  the  noise 
measurements  performed  using  a  hydrophone.  Figure  8.16  shows 
the  peak  pressures  plotted  against  6  for  each  of  the  bubble 
configurations.  One  can  see  that  the  highest  pressures  were 
recorded  within  a  small  range  of  8,  (0.7  <  8  <  0.13).  Also,  the 
peak  pressure  increases  with  higher  bubble  density  configuration 
(higher  N  here)  as  the  theory  predicts.  This  data  again 
suggests  the  existence  of  an  optimum  value  for  8  fairly  close  to 
that  mentioned  above.  At  this  optimum  value,  the  effect  of  high 
field  pressures  and  large  jet  penetration  seem  to  exist 
simultaneously . 

For  all  of  the  tests  considered  above,  the  designated 
original  center  of  the  bubble  was  seen  to  move  more  or  less 
significantly  toward  the  center  of  the  cloud.  This  is  shown 
graphically  in  Figures  8.17-8.20,  which  show  the  motion  of 
the  bubble  center.  The  distance,  6  =  R3  -  R2,  is  shown  versus 
time  for  the  different  bubble  configurations.  One  can  see  that 
very  little  bubble  "center"  motion  was  occuring  before  the  first 
collapse.  However,  once  collapsed,  the  bubbles  both  grow  again 
and  move  rapidly  toward  the  center  of  the  cloud.  The  shift 
toward  the  center  is  greater  for  the  higher  density  bubble 
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configuration.  This,  again,  relates  directly  to  the  results 
discussed  above,  since  large  values  for  5  correspond  to  smaller 
values  for  R2,  or  larger  penetratio  t  of  the  re-entering  jets. 
Figures  8.21  and  8.22  show  typical  sequences  of  pictures  taken 
from  high  speed  movies.  Bubble  deformation  and  motion  towards 
the  cone  apex  can  be  clearly  observed. 

8.3  Conclusions 


A  very  extensive  set  of  experiments  was  run  to  study  the 
influence  of  bubble  number  and  spacing  on  the  dynamics  of  a 
simulated  multibubble  system.  High  speed  photography  and  noise 
measurements  allowed  us  to  follow  with  precision  the  dynamics  of 
spark-generated  bubbles.  However,  as  already  discussed  in 
Section  4.4,  a  major  impediment  to  confident  conclusions  is  the 
repeatibil ity  of  the  tests.  Indeed,  the  greatest  source  of 
error  in  the  interpretation  of  the  results  is  the  amount  cf  gas 
in  the  bubble  which  is  not  easily  controlled.  With  this  in 
mind,  there  is  a  strong  indication  that  the  collapse  of  the 
simulated  cloud  is  very  intense  for  an  optimum  value  of  e,  the 
ratio  of  characteristic  bubble  size  to  bubble  interdistance. 

When  e  increases  (starting  from  zero)  the  bubble  period  of 
oscillation  increases  first,  and  the  collapse  pressure  decreases 
slightly.  Later,  a  strengthening  of  the  implosion  occurs 
accompanying  an  increase  of  bubble  collapse  and  a  shortening 
of  the  bubble  period.  These  results  match,  qualitatively  at 
least,  the  predicted  theoretical  behavior.  As  e  keeps 
increasing,  the  experiments  show  other  pressure  peaks,  or  a 
collapse  period  minimum  for  certain  particular  values  of  e. 

This  result,  at  first  unexpected,  seems  to  correspond  to 
theoretical  predictions  for  the  case  of  high  gas  content  in  the 
bubbles  (see  Section  4.4). 


9.0  ACOUSTIC  DAMPING  OF  A  PRESSURE  FIELD  BY  A  BUBBLE  SCREEN 


9 . 1  Objectives 

In  order  to  complement  the  spark-generated  bubble 
experience  where  simulated  symmetrical  cloud  configurations  were 
submitted  to  strong  oscillations,  the  acoustical  behavior  of  a 
bubble  cloud  (screen)  was  investigated.  Here  no  symmetry  of  the 
bubble  configuration  or  size  was  involved  and  the  imposed 
pressure  oscillations  were  moderate.  The  main  objective  of  this 
preliminary  study  was  to  find  correlations  between  bubble  size 
distribution  and  bubble  concentration  and  the  transmission  of 
various  frequency  sound  waves  through  the  cloud.  From  single 
bubble  studies  it  is  known  that  there  is  a  linear  relation 
between  bubble  size  and  bubble  natural  frequency.  The  acoustic 
behavior  of  the  bubble  is  also  known  to  be  a  function  of  the 
relative  magnitudes  of  the  natural  frequency  and  the  frequency 
of  the  excitation.  Acoustic  resonance  occurs  when  both 
frequencies  are  the  same.  In  the  present  tests  interaction 
between  bubbles  of  different  sizes  in  the  screen  come  into  play, 
and  the  overall  results  should  account  for  these  effects. 

9 . 2  Description  of  Experiment 

To  experimentally  determine  a  relationship  between  bubble 
size,  density,  and  acoustic  damping,  a  number  of  tests  were 
conducted  varying  the  density  of  air  bubbles  in  water  and 
recording  the  change  in  amplitude  of  acoustic  signals  across  a 
path  of  the  bubbles.  The  signals  ranged  in  frequency  from  5kHz 
to  70kHz.  Photographs  were  taken  at  the  time  of  each  test  to 
determine  numerical  values  for  bubble  size  distribution. 
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The  tests  were  conducted  in  a  large  plexiglass  tank 
approximately  5  feet  x  4  feet  x  1  foot.  The  tank  was  filled 
with  water  to  a  depth  of  3.5  feet.  Four  porous  nylon  tubes, 
HYDROPERM  tubes,  5/8  inch  diameter,  11  inches  long,  and  2 
inches  apart  were  placed  in  the  bottom  of  the  tank  to  generate 
bubbles.  The  approximate  pore  size  for  these  tubes  was  1-3 
microns.  All  but  the  top  portion  of  the  circumference  of  the 
nylon  tubes  was  sealed  off  in  order  to  avoid  the  coalescence  of 
bubbles  from  the  sides  and  bottom  of  the  tubes.  Each  porous 
tube  was  connected  to  a  common  pipe  at  one  end.  This  pipe  was 
then  attached  to  a  compressed  air  line  to  regulate  and  measure 
the  pressure  of  the  air  being  forced  into  the  porous  tubes.  For 
the  purpose  of  these  tests,  the  air  pressure  in  the  tubes  ranged 
from  6  psig  to  16  psig. 

Two  hydrophones  were  placed  in  the  water  9  inches  apart  and 
21  inches  below  the  surface  of  the  water.  One  hydrophone  served 
as  a  signal  transmitter,  and  the  other  as  a  signal  receiver. 

The  transmitting  hydrophone  was  a  rubber  booted,  ceramic 
cylinder  type  with  a  practically  flat  response  up  to  15kHz  and 
±3dB  deviation  between  15  and  30  kHz.  A  KSP  UT-114  hydrophone, 
having  a  flat  response  up  to  30  kHz  and  ±3dB  oscillations  up  to 
50kHz,  was  used  as  a  signal  receiver.  Signals  were  produced  at 
various  frequencies  using  a  Wavetek  Model  113  frequency 
generator,  amplified  with  a  McIntosh  MC75  audio  amplifier,  then 
transmitted  by  the  hydrophone.  The  signal  was  then  picked  up  by 
the  receiving  hydrophone,  where  it  was  first  amplified  then 
filtered  to  eliminate  excessive  low  frequency  noise  (f<  500 
Hz).  The  received  signals  were  monitored  with  a  Hewlett  Packard 
3580A  frequency  analyzer,  (maximum  frequency  of  50  kHz),  a 
Keithley  197  digital  multimeter,  and  a  Nicolet  2090  digital 
oscilloscope.  Figure  9.1  is  a  diagram  of  this  experimental 
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set-up.  Figure  9.2  shows  a  photograph  of  the  actual  tank  itself 
with  the  production  of  bubbles  by  16  psig  air  pressure  being  fed 
through  the  porous  tubes.  It  shows  the  arrangement  and  location 
of  the  two  hydrophones  with  bubbles  passing  between  them. 

Signals  of  various  frequencies  were  first  generated  with  no 
air  bubbles  being  produced  in  the  tank.  The  strengths  of  these 
signals  as  measured  by  the  multimeter  were  recorded  for  later 
comparison  with  those  transmitted  across  a  stream  of  air 
bubbles.  In  this  manner,  the  amount  of  signal  damping  which 
occurred  at  each  frequency  due  to  the  bubbles  could  be 
determined.  This  could  then  be  related  to  a  characteristic 
bubble  size  distribution  by  analyzing  the  photographs  taken  at 
the  time  of  each  test. 

9 . 3  Experimental  Results  and  Discussion 

Figure  9.3  shows  the  results  of  some  of  the  early  tests  run 
with  the  above  described  apparatus.  The  values  for  the  air 
pressure  ranged  from  8  psig  to  16  psig  and  the  results  were 
somewhat  encouraging.  The  highest  damping  (ratio  of  rms  signal 
emitted  minus  signal  detected  to  signal  emitted  )  occurred  at 
the  higher  bubble  densities,  as  is  expected.  However,  the  sharp 
drop  in  the  amount  of  damping  which  occurred  at  around  40  kHz 
for  most  of  the  tests  could  not  be  explained  from  the 
characteristics  of  the  bubbles  themselves. 

After  some  searching,  it  was  discovered  that  the  behavior 
of  the  curves  shown  in  Figure  9.3  was  due,  in  fact,  more  to  the 
response  of  the  hydrophones  than  to  the  size  of  the  bubbles.  To 
correct  this,  the  receiving  hydrophone  was  replaced  with  the 
one  described  above;  whereas  before,  this  hydrophone  had  been 
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one  similar  to  the  transmitting  hydrophone.  With  two  practically 
identical  hydrophones  being  used,  the  variations  of  the 
frequency  response  which  were  occurring  at  the  higher 
frequencies  were  being  overly  amplified.  Figure  9.4  shows  the 
calibrated  frequency  responses  of  both  hydrophones. 

With  the  new  receiving  hydrophone  in  place,  a  new  series  of 
tests  was  conducted  to  determine  the  range  of  frequencies  at 
which  the  most  acoustic  damping  would  occur.  Figure  9.5  shows 
the  results  of  a  test  run  at  an  air  pressure  of  10  psig.  The 
plot  presents  both  the  percentage  damping  which  occurred  with 
bubbles  and  the  amplitude  of  the  received  signal  with  no  bubbles 
being  generated  in  the  tank.  It  suggests  that  most  of  the 
damping  is  occurring  with  frequencies  in  the  range  of  20-50kHz. 
But  the  response  of  the  signals  with  no  bubbles  is  still 
suspicious  with  sharp  spikes  at  about  35  and  50kHz.  To 
determine  any  effect  this  may  still  be  having  on  the  shape  of 
the  damped  signal,  a  voltmeter  was  added  to  the  experimental 
set-up  to  monitor  the  input  signal  to  the  transmitting 
hydrophone.  A  new  test  was  run  maintaining  the  output  signal 
constant  in  the  absence  of  bubbles  and  recording  the  input 
voltage  necessary  to  do  so  over  the  range  of  frequencies  of 
interest.  Then,  with  bubbles  being  produced,  the  measured  input 
voltage  at  each  particular  Crequency  was  matched  to  that 
recorded  with  no  bubbles.  The  corresponding  amplitude  of  the 
output  signal  was  then  recorded.  The  results  of  this  test  are 
shown  in  Figure  9.6.  One  can  see  that,  although  the  maximum 
amount  of  damping  is  slightly  lower  than  that  shown  in  Figure 
9.7,  the  overall  behavior  is  still  very  much  the  same. 

Figures  9.7  and  9.8  are  photographs  taken  from  the  front  of 
the  tank  of  the  bubbles  at  10  and  16  psig  air  pressures, 
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respectively .  The  photos  suggest  that  there  exists  a  range  of 
bubble  sizes  and  that  the  sizes  changed  very  little  with  the 
increasing  air  pressure  in  the  tubes.  This  would  explain 
qualitatively  the  results  presented  above.  The  distribution  of 
bubble  sizes  is  responsible  for  the  range  of  frequencies  over 
which  acoustical  damping  of  the  signals  has  occurred.  Similar 
results  for  tests  of  various  air  pressures  in  the  porous  tubes 
suggest  that  the  bubble  sizes  were  similar  for  all  the  tests,  as 
is  evidenced  by  these  two  photographs. 

In  the  two  photographs,  Figures  9.7  and  9.8,  can  be  seen 
three  small  wires,  5,  10,  and  15  mils,  used  for  scaling  purposes 
when  analyzing  the  photographs  for  bubble  size  distribution. 
Figure  9.9  shows  the  spectrum  of  bubble  radii  as  was  determined 
from  the  photo  in  Figure  9.7.  It  shows  that  most  of  the  bubble 
radii  fall  within  a  range  of  0.05  mm  to  0.30  mm.  This  is  only  a 
rough  estimate  of  the  sizes  due  to  the  fact  that  the  bubbles  are 
too  small  to  be  accurately  measured.  However,  when  comparing  to 
theoretical  considerations,  these  results  seem  to  be  accurate. 
Figure  9.10  shows  graphically  the  theoretical  relationship 
between  bubble  size  and  resonance  frequency  for  an  isothermal 
compression  law.  From  it,  one  can  see  that  the  particular 
bubble  size  range  given  above  corresponds  theoretically  to 
frequencies  of  10-55kHz.  This  is  a  very  good  match  with  the 
experimental  results  presented  above.  The  wide  range  of 
frequencies  that  are  significantly  damped  ( 1 5  <  f  <  60  KHz,  in 
Figure  9.6)  corresponds  to  the  measured  bubble  sizes.  However, 
as  seen  in  Figure  9.9,  the  acoustical  measurements  seem  to 
indicate  a  higher  content  of  smaller  bubbles  than  visually 
observed.  This  could  be  due  to  higher  modes  of  bubble 
oscillations  or  to  deviations  from  the  simple  resonance 
frequency-bubble  size  relationship  shown  in  Figure  9.10. 
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To  further  verify  the  relationship  between  bubble  size  and 
acoustically  attenuated  frequencies,  the  experimental  apparatus 
described  earlier  was  slightly  altered  in  two  separate  ways  in 
an  effort  to  increase  the  size  of  the  bubbles  being  produced  in 
the  tank.  In  addition  to  these  alterations,  bubbles  were  also 
produced  by  electrolysis  with  electrical  wires.  These  bubbles 
proved  to  be  too  small  for  our  purposes  (i.e.,  natural  frequency 
too  high  for  the  instrumentation).  The  first  change  consisted 
of  placing  in  the  tank  a  set  of  stainless  steel  sintered  tubes, 
(pore  size  =  5  microns),  in  much  the  same  arrangement  as  the 
original  nylon  tubes.  These  tubes  were  again  sealed  along  all 
but  the  very  top  portion  of  the  circumference.  Signals  of 
frequencies  ranging  from  1-50  kHz  were  generated  in  the  tank  and 
the  results  showed  that  a  considerable  amount  of  damping  was 
occurring  at  practically  all  the  frequencies  within  this 
spectrum.  While  these  results  were  useful  in  showing 
attenuation  of  lower  frequencies  than  the  HYDROPERM  tubes,  they 
did  not  exhibit  the  definite  shift  in  frequency  bandwidth  as  had 
been  hoped  for  with  this  new  apparatus.  Apparently,  with  the 
sintered  tubes,  bubbles  with  a  very  wide  range  of  sizes  were 
being  produced. 

As  a  second  attempt  to  produce  bubbles  of  a  noticeable 
difference  in  size  from  those  produced  by  the  nylon  tubes,  an 
aluminum  tube  with  drilled  holes,  0.006  inch  diameter,  was  used 
to  generate  bubbles  at  the  bottom  of  the  test  tank.  Figures 
9.11  and  9.12  show  photographs  of  the  bubbles  produced  from  this 
apparatus  taken  at  air  pressures  of  3  and  7  psig,  respectively. 
Again,  inspection  of  these  two  figures  shows  that  the  size 
distribution  of  the  bubbles  being  produced  at  these  pressures 
changed  very  little  if  at  all.  As  before,  a  size  distribution 
was  determined  from  the  photographs  and  is  given  in  Figure 
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9.13.  This  figure  suggests  that  most  of  the  bubble  radii  were 
within  a  range  of  1.22-1.45  mm.  From  Figure  9.10,  the  damped 
frequencies  for  this  range  of  bubble  sizes  should  fall  between 
1.9  and  2.2  kHz.  Figure  9.14  shows  the  actual  experimental 
results  of  the  tests  run  with  these  bubbles.  Problems  related  to 
data  reproduction  at  various  dates  can  be  observed  in  this 
figure.  Peak  damping  occurred  at  about  2  kHz  and  a  much 
narrower  bandwidth  of  frequencies  was  being  attenuated  by  the 
presence  of  the  bubbles.  This  corresponds  very  well  with 
theoretical  considerations  since  a  spectrum  of  larger  bubble 
radii  (>1  mm)  results  in  very  little  variation  in  resonant 
frequency,  whereas  a  range  of  bubble  radii  <0.5  mm  will  produce 
a  much  larger  bandwidth  of  resonant  frequencies.  For  the  purpose 
of  these  tests,  no  results  could  be  obtained  at  frequencies 
lower  than  2  kHz  because  of  an  increased  noise  to  signal  ratio 
at  these  frequencies,  but  the  results  do  show  the  shift  toward 
lower  frequencies  with  increased  bubble  size  as  was  expected. 
This  set  of  tests  give  some  preliminary  but  very  useful 
information  on  the  study  of  acoustic  wave  propagation  through  a 
bubbly  medium  . 
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10.0  CONCLUSIONS  AND  SUMMARY 

The  main  achievement  of  the  study  described  in  this  report 
is  the  extension  of  the  study  of  single  bubble  dynamics  to 
multibubble  dynamics.  Within  the  restriction  that 
characteristic  bubble  size  is  small  compared  to  interbubble 
distance,  the  dynamics  of  a  bubble  cloud  was  investigated. 
Dynamic  effects  in  the  absence  of  heat  and  mass  transfer  were 
first  considered.  The  matched  asymptotic  expansion  method  was 
used  and  enabled  derivation  of  differential  equations  for  the 
bubble  shape  components  at  the  various  orders  of  approximation. 
These  were  solved  numerically  using  a  multi-Runge-Kutta 
procedure,  and  the  corresponding  numerical  code  was  implemented. 

Later,  heat  then  mass  transfer  were  accounted  for.  When, 
during  a  significant  portion  of  the  imposed  pressure  field 
history,  the  liquid  is  superheated  (its  temperature  is  higher 
than  the  liquid  vaporization  temperature  at  the  imposed 
pressure)  or  gas  supersaturated  (the  concentration  of  dissolved 
gas  is  higher  than  the  saturation  concentration  at  the 
noncondensible  gas  pressure  in  the  bubbles),  heat  and  mass 
transfer  at  the  interface  occur  at  a  rate  high  enough  to 
interfere  with  the  dynamics  of  the  bubble  growth  or 
oscillations.  In  that  case,  the  dynamic  equation  of  a  bubble  in 
the  cloud  is  coupled  to  the  heat  or  mass  transfer  equation  by 
the  value  of  one  (or  both)  of  the  partial  pressures  of  the 
liquid  vapor  and  noncondensible  gas  inside  the  bubble  which  are 
then  transfer  dependent.  We  have  solved  this  problem  when  (as 
for  single  isolated  bubbles)  the  liquid  properties  or  conditions 
are  such  that  the  variation  of  the  liquid  temperature  or  the  gas 
concentration  in  the  liquid  occurs  primarily  in  a  thin  boundary 
layer  at  the  bubble  wall.  In  that  case  an  integral  equations 
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relates  bubble  dynamics  and  heat  or  mass  transfer  at  the  bubble 
interface.  Two  numerical  codes  were  developed  to  deal  with  the 
two  transfer  problems. 

Other  aspects  of  the  problem  were  also  considered  in  this 
study  but  were  not  pursued  as  deeply  as  the  multibubble  approach 
due  either  to  limited  time  and  resources  or  to  their  showing 
little  promise  based  on  initial  feasibility  efforts.  The  most 
promising  of  these  approaches  have  been  summarized  in  this 
report.  Some  others,  such  as  the  determination  of  the  unsteady 
pressure  field  generated  behind  an  oscillating  sheet  cavity  have 
not  been  addressed  in  this  report.  The  task  turned  out  to  be 
much  more  involved  than  originally  thought  if  a  significant  new 
contribution  was  to  be  achieved.  The  effort  was  redirected 
towards  a  more  thorough  analysis  of  the  other  tasks. 

A  continuum  medium  approach  was  developed  for  the  general 
description  of  the  two  phase  medium  composing  the  cloud  that 
shows  promise.  As  do  most  continuum  approaches,  this  method 
applies  to  the  case  of  low  void  fractions.  It  enables 
expression  of  the  conservation  equations  in  the  cloud  medium. 

The  micromorphic  method  used  consists  in  decomposing  the  medium 
into  cells  in  which  the  bubble  dynamics  are  considered.  One 
shortcoming  is  the  appearance  of  an  additional  term  arising  from 
the  averaging  procedure  which  depends  on  the  geometry  of  the 
cell.  Dependence  of  the  equations  on  an  inverse  power  of  the 
void  fraction  indicates  they  are  invalid  for  extremely  low  void 
fractions.  These  difficulties  are  classic  with  averaging 
methods  and  can  now  be  dealt  with  using  newly  developed 
renormalization  techniques.  This  should  be  the  next  step 
towards  the  improvement  of  this  continuum  theory. 
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A  fundamental  experimental  program  was  conducted  in 
parallel  with  the  theoretical  effort  described  here.  Systematic 
observations  of  bubble  dynamics  using  high  speed  photography  and 
measurements  of  acoustical  pressures  were  conducted. 
Spark-generated  bubbles  were  produced  at  precise  locations  near 
complex  solid  wall  geometries.  By  application  of  the  method  of 
images,  the  bubble  behavior  near  the  solid  boundaries  is  similar 
to  that  in  the  presence  of  a  set  of  bubble  images  in  the  wall. 

By  varying  the  wall  geometry  the  number  of  simulated  interacting 
bubbles  can  be  varied.  The  interbubble  distance  was  set  by 
controlling  the  distance  between  the  electrode  tips  and  the 
solid  walls.  The  size  of  the  bubble  generated  was  varied  by 
changing  the  discharge  voltage  and  the  tank  ambient  pressure. 
Observations  confirmed  and  complemented  the  theory  since  cases 
of  both  low  and  high  void  fractions  could  be  studied.  Within 
the  experimental  error,  (mainly  due  to  a  lack  of  control  of  the 
amount  of  gas  in  the  bubble),  preferred  void  fractions  (or 
ratios  of  bubble  characteristic  radius  to  interbubble  distance) 
which  generated  the  strongest  bubble  collapse  and  the  highest 
pressures  were  observed.  These  results  qualitatively  confirmed 
the  theoretical  results. 

A  set  of  experiments  of  a  more  preliminary  nature  concerned 
sound  wave  propagation  through  a  sheet  of  bubbles.  Bubble 
sheets  were  generated  in  a  large  tank  using  microporous, 
sintered,  and  drilled  tubes.  Sound  pressure  waves  were  then 
emitted  on  one  side  of  the  sheet  using  an  active  hydrophone. 

The  transmitted  wave  was  detected  with  another  hydrophone.  By 
varying  the  frequency  of  the  emitted  signals  and  the 
characteristics  of  the  bubble  generator,  correlations  could  be 
made  between  the  attenuation  of  the  acoustic  wave  through  the 
bubble  screen  and  the  size  distribution  of  the  bubbles. 
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Several  numerical  codes  were  implemented  as  a  result  of  the 
multibubble  cloud  theory  developed  in  this  program.  These 
programs  consider  several  aspects  of  the  problem.  CLDMAIN 
neglects  heat  and  mass  transfer  and  allows  the  user  to  study  the 
behavior  of  a  symmetrical  bubble  cloud  configuration.  In  its 
latest  verison,  CLDMAIN2,  several  pressure  field  shapes  can  be 
investigated  —  e.g.,  sudden  pressure  rise,  pressure  drop  for  a 
finite  time  period,  pressure  field  due  to  a  venturi  or  over  a 
foil  shape,  ...  etc.  The  error  at  each  step  of  computation  can 
be  bounded  by  an  imposed  value.  This  procedure  is  used  to 
continuously  adjust  the  calculation  time  step  to  minimize  both 
errors  and  total  computation  time.  This  code  has  been 
extensively  exercised  to  investigate  bubble  cloud  behavior.  Due 
to  the  symmetry  of  the  cloud,  very  dramatic  results  can  be 
observed.  A  cumulative  effect  is  observed  which  tends  to 
reinforce  any  pressure  deviations  from  the  ambient  value  due  to 
bubble  dynamics.  For  instance,  during  the  initial  bubble  growth 
period  the  positive  deviations  in  the  presence  of  several 
bubbles  reinforce  each  other  to  slow  down  the  initial  bubble 
growth.  This  effect  can,  however,  be  reversed  later  when  the 
deviation  pressure  becomes  negative  and  can  lead  to  large  bubble 
sizes.  The  most  dramatic  effect  is  observed  during  bubble 
collapse.  It  is  known  that  large  pressure  rises  are  generated 
at  the  end  of  the  collapse.  These  pressures  are  increased  when 
several  bubbles  are  present  and  also  reinforce  the  collapse  of 
each  individual  bubble.  Therefore,  collective  bubble  collapse 
can  generate  pressures  orders  of  magnitude  higher  than  those 
produced  by  single  bubble  collapse.  The  cumulative  effect  is 
mainly  due  to  the  fact  that  each  bubble  ends  its  collapse  under 
the  influence  of  pressures  generated  by  the  collapse  of  the 
other  bubbles  orders  of  magnitude  higher  than  the  imposed 
ambient  pressure.  This  would  tend  to  explain  the  observed  high 
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erosion  intensities  and  the  bending  of  some  propeller  trailing 
edges . 

When  the  bubble  cloud  geometry  is  asymmetrical  the 
pressures  generated  by  the  various  bubbles  are  not  necessarily 
in  phase.  Cumulative  effects  are  therefore  reduced  in  intensity 
and  the  very  large  increases  in  the  collapse  energy  found  with 
the  symmetrical  CLDMAIN  code  are  tempered.  A  similar  result  is 
obtained  when  bubble  initial  sizes  are  not  identical  or  when  a 
finite  sound  speed  is  taken  into  account  to  delay  pressure 
propagation  from  one  bubble  to  another.  Two  main  codes  have 
been  developed  to  deal  with  these  asymmetries:  MULTIBBL  and 
MULTICOMP.  For  both  codes  the  user  inputs  bubble  size  and 
distribution,  and  the  program  computes  the  individual  bubble 
dynamics  under  the  influence  of  the  other  bubbles  in  the  cloud. 
In  addition,  MULTICOMP  accounts  for  time  delays  between  the 
arrival  of  pressure  signals  from  different  bubbles  based  on  a 
finite  sound  speed  imposed  by  the  user. 

For  heat  and  mass  transfer,  two  separate  codes  have  been 
implemented:  HOTCLDS  and  GASCLD.  Both  consider  a  symmetrical 

cloud  configuration  and  couple  the  multibubble  dynamics  problem 
with  the  problem  of  heat  or  mass  transfer  at  the  bubble  liquid 
interface.  Mass  and  heat  transfer  are  both  seen  to  reduce 
bubble  growth  rate.  Bubble  interactions  also  play  a  significant 
role  in  bubble  growth  —  mainly  by  modification  of  the  pressure 
field.  This  interaction  results  in  a  reduced  growth  rate  which 
makes  the  bubble  radius  at  any  given  time  smaller  than  would  be 
found  for  an  isolated  bubble.  The  temperature  drop  at  the 
bubble  wall  in  a  multibubble  configuration  is  smaller  than  would 
be  found  for  an  isolated  bubble.  The  dissolved  gas  concentra¬ 
tion  in  the  liquid  exhibits  similar  behavior.  The  effective  gas 
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compression  law  is  seen  then  to  vary  rapidly  from  an  adiabatic 
behavior  at  the  beginning  of  the  growth  to  an  isothermal 
behavior. 

In  conclusion,  we  have  conducted  an  in-depth  study  on  the 
fundamental  aspects  of  multibubble  cloud  cavitation.  The 
results  of  this  study  explain  the  intense  erosion  associated 
with  the  phenomena.  Several  codes  were  developed  and 
implemented  as  a  result  of  the  theoretical  study  and  can  be 
used  as  tools  for  further  study  of  various  aspects  of  the 
problem  under  particular  conditions  of  interest. 
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TABLE  1 .  VALUES  OF  THE  NUMERICAL  CONSTANTS 


USED  IN  THE  COMPUTATIONS 
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FIGURE  1.1  -  TYPES  OF  CAVITATION  ON  FOIL 
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FIGURE  3.1  -  MULTIBUBBLE  INTERACTION  EQUIVALENCE  CONCEPT 
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CIGURE  4.1  -  BUBBLE  CLOUD  DYNAMICS  PROCRAM  FLOW  CHART. 

EQUI-SIZED  SYMMETRICALLY  DISTRIBUTED  BUBBLES 
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FIGURE  4.4  -  PRESSURES  GENERATED  AT  A  POINT  M  BY  TWO  COLLAPSING 
BUBBLES:  BUBBLE  A,  SMALLER  OR  FARTHER  AWAY;  BUBBLE 
B,  LARGER  OR  CLOSER 


FIGURE  4.8  -  INFLUENCE  OF  THE  PRESSURE  DROP  DURATION  ON  THE 
BUBBLE  WALL  MOTION  TOWARD  THE  CLOUD  CENTER, 


(Tilt 


DISTANCE  r  FROM 


BUBBLES  E  P  S  =  0 . 1 


FIGURE  4.13  -  INFLUENCE  OF  THE  INITIAL  GAS  PRESSURE  ON  THE  COLLAPSE 

OF  A  FOUR-BUBBLE  CLOUD.  C,  =  2.707,  C,  =  1.914,  C„  =  0.85 


FIGURE  4.14  -  INFLUENCE  OF  e  ON  THE  COLLAPSE  OF  A  FOUR-BUBBLE  CLOUD 
C,  =  2.707,  C„  =  1.914,  C,  =  0.854,  P  =  0.1 


NONDIMENS  I  ON  AL  RADIUS,  R  Ir 
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NORMALIZED  TIME,  t/RQ  V  p 7  A  P 


FIGURE  4.16  -  INFLUENCE  OF  e  ON  THE  COLLAPSE  OF  A  TWELVE-BUBBLE 

CLOUD,  C1  =  8.616,  C2  =  4.53,  C3  -  0.41,  P  =  0.10 

90 
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FIGURE  4.17  -  INFLUENCE  OF  e  ON  THE  TIME  TO  FIRST  COLLAPSE  FOR  A 

FOUR-BUBBLE  CONFIGURATION  AND  VARIOUS  INITIAL  GAS 
PRESSURES 


FIGURE  4.18  -  INFLUENCE  OF  e  ON  THE  TIME  TO  FIRST  COLLAPSE  FOR  AN 

EIGHT-BUBBLE  CONFIGURATION.  P  =  0.05 
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FIGURE  4.19  -  INFLUENCE  OF  £  ON  THE  TIME  TO  FIRST  COLLAPSE  FOR  A 

TWELVE-BUBBLE  CONFIGURATION  FOR  VARIOUS  INITIAL  CAS 
PRESSURES 


FIGURE  4.21  -  CRITICAL  VELOCITY  AT  VENTURI  THROAT  AS  A  FUNCTION  OF  THE 
INITIAL  BUBBLE  RADIUS.  STATIC  EQUILIBRIUM  CASE,  K  =  1.4 


POLYTROPIC  CONSTANT 


BUBBLE  RADIUS, 
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FIGURE  4.25 


THROAT  VELOCITY,  V  ,  m/s 

MAXIMUM  BUBBLE  RADIUS  AS  A  FUNCTION  OF  THE 
THROAT  VELOCITY.  INFLUENCE  OF  BUBBLE  SIZE 
AND  NUMBER 


FIGURE  4.27  -  BUBBLE  GROWTH  AND  COLLAPSE  IN  THE  VENTURI 
SYSTEM  AT  THREE  DIFFERENT  INITIAL  BUBBLE 


VENTURI  BUBBLE  DYNAMICS 


FIGURE  4.28  -  COMPARISON  BETWEEN  DYNAMICS  AND  STATIC  CRITICAL  VELOCITY 
AT  VENTURI  THROAT  AS  A  FUNCTION  OF  INITIAL  BUBBLE  RADIUS 


EPS=0 . 1 ,  WEB= 1 00 . ,  PGZB=0 . 2 ,  Lo=10,  Rbo= 1 
R 1 =R2= 1 ,  R3=0 . 3 ,  MAJOR  RADIUS-AR 


FIGURE  4.29  -  COLLAPSE  HISTORY  OF  A  THREE  BUBBLE  CONFIGURATI 
RESPECTIVE  INITIAL  NONDIMENS IONAL  RADII  OF  1,  1  A 
COMPARISON  WITH  THE  ISOLATED  BUBBLE  CASES 


EPS=0 . 1 ,  WEB=100..  PGZB=0 . 2 ,  Lo=10,  Rbo= f 
R1 -  1  R2=R3=0 . 3 ,  MAJOR  RADIUS- AR 


FIGURE  4.30  -  COLLAPSE  HISTORY  OF  A  THREE  BUBBLE  CONFIGURATION. 

RESPECTIVE  INITIAL  NON  DIMENSIONAL  RADII  OF  1,  0.3 
AND  0.3  COMPARISON  WITH  THE  ISOLATED  BUBBLE  CASES 


.  WEB= 1 00 . ,  PCZB=0.2, 

BUBBLE  SYSTEHS  COMPARSION,  MAJOR 
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FIGURE  4.31  -  COMPARISON  OF  THE  RADIUS  HISTORY  FOR  A  SYMMETRICAL  12-BUBBLE 
CONFIGURATION  (BUBBLES  CENTERED  ON  SPHERE  R  =  15),  AND  A  2- 
LAYER  6-BUBBLE  CONFIGURATION  (BUBBLES  CENTERED  ON  SPHERES 
R  =  7.5  AND  R  =  15).  BOTH  CONFIGURATIONS  HAVE  THE  SAME  GLOBAL 
VOID  FRACTION 


FIGURE  4.33  -  BUBBLE  SIZES  HISTORY  FOR  A  RANDOMLY  SIZED  AND  RANDOMLY 

DISTRIBUTED  6-BUBBLE  CONFIGURATION  (SEE  TABLE  FIGURE  4.34) 
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BUBBLE  NO. 

COORDINATES 

RADIUS 
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FIGURE  4.34 


RANDOM  BUBBLE  DISTRIBUTION  AND  SIZES 
(Corresponds  to  results  on  Figure  4.33) 


BllBBLES-6  BUBBLES  WITH  WALL  DISTANCED.  MINOR  RADIUS-TOP 
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FIGURE  4.35  -  GROWTH  AND  COLLAPSE  OF  A  6-BUBBLE  CLOUD  NEAR  A  SOLID 
WALL.  RADIUS  HISTORY  FOR  THE  VARIOUS  BUBBLES:  NO.  1. 
CLOSEST  TO  THE  WALL.  NO.  3.  FARTHEST  FROM  THE  WALL 
NO.  6.  FOUR  BUBBLES  IN  A  PLANE  PARALLEL  TO  THE  WALL. 
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FIGURE  5.2  -  RATIO  OF  NON-SPHERICAL  TO  SPHERICAL  COMPONENT  OF  THE 

RADIUS  VERSUS  TIME  FOR  A  BUBBLE  GROWING  IN  AN  N-BUBBLE 
CLOUD.  W  =  100.  R  =  0.0001  m,  e  =  0.1 


DATA  INPUT  PROCRAM,  HOT  IN.  I - H  DATA  READING 


FIGURE  5.3  -  PROGRAM  FLOW  CHART  OF  HOTCLOS.  SYMMETRICAL  CLOUD  DYNAMICS 

ACCOUNTING  FOR  HEAT  LXCHAN< .L  AND  PHASE  CHANGE 
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FIGURE  5.6  -  VARIATION  WITH  TIME  OF  BUBBLE  SHAPE  CHARACTERISTICS 
AT  TWO  SELECTED  TIMES,  Po  =  2  atm,  Pjnf  =  0.5  atm,  e  =  0.05 

Rfa  =  0.01  m,  RQ  =  2.5  x  10_5m,  Tfa  =  809. 3°C,  T,nf  =  1177. 71°C 
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FIGURE  5.10  -  INFLUENCE  OF  IMPOSED  PRESSURE  P.  ,  ON  THE  BUBBLE-WALL  TEMPERATURE 
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FIGURE  5.13  -  INFLUENCE  OF  INITIAL  BUBBLE  SIZE  ON  BUBBLE-WALL-TEMPERATURE 

HISTORY.  PQ  =  2  atm,  Pjnf  =  1  atm.  e  =  0.2,  Rfa  =  0.01  m,  Tfa  =  879. 9°C: 

A,  R  =  2.5E  -  6;  B,  2.5E  -  5;  C,  1.0E  -  4,  D,  °1.0E  -  3 
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FIGURE  6.1  -  PROGRAM  FLOW  CHART  FOR  GASCLO.  SYMMETRICAL  CLOUD  DYNAMICS 

ACCOUNTING  FOR  NON-CONDENSIBLE  CAS  TRANSFER  ACROSS  THE 
BUBBLE  WALL 


Vv 


FIGURE  6.6  -  VARIATION  OF  EFFECTIVE  GAS  POLYTROPIC  COEFFICIENT 

DURING  BUBBLE  GROWTH  FOLLOWING  A  SUDDEN  DEPRES¬ 
SURIZATION  ( P  =5  atm,  P.  ,  =  0.01  atm,  R.  =  1  mm. 


FIGURE  6.7  -  INFLUENCE  OF  INITIAL  DISOLVED  GAS  CONCENTRATION 

ON  EFFECTIVE  GAS  POLYTROPIC  COEFFICIENT  (P  =  3.7 


FIGURE  6.8  -  INFLUENCE  OF  INITIAL  DISSOLVED  GAS  CONCENTRATION  ON 
EFFECTIVE  GAS  POLYTROPIC  COEFFICIENTS  (P„  =  5  atm. 
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FIGURE  8.3  -  PRESSURE  VARIATION  VS.  TIME  FOR  SPARK  INDUCED 
BUBBLE  2i  INCH  FROM  CONE  TIP 


FIGURE  8.4  -  FREQUENCY  TRANSFORM  OF  DATA  SHOWN  IN  FIGURE  26 
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B  :  INITIAL  BUBBLE  CENTER; 
OR  ELECTRODE  CAP 


B'  :  NEW  POSITION  OF  BUBBLE 
CENTER 

O  :  BUBBLE  CLOUD  CENTER; 

OR  PYRAMID  APEX 

FIGURE  8.5  -  PRINCIPAL  GEOMETRICAL  CHARACTER¬ 
ISTICS  OF  A  BUBBLE  DURING  ITS  MOTION 


FIGURE  8.7  -  VARIATIONS  OF  R2  VERSUS  TIME  FOR  DIFFERENT  VALUES 
OF  8  .  TESTS  CONDUCTED  WITH  3-DIMENSIONAL  CORNER 


PYRAMID 


TESTS  CONDUCTED  WITH  PENTAGONAL  PYRAMID, 


VS.  BETA 


THE  VARIOUS  BUBBLE  CONFIGURATIONS 


TESTS  CONDUCTED  WITH  3-DIMENSIONAL  CORNER 


FIGURE  8.13  -  VARIATION  OF  R3  VERSUS  TIME  FOR  DIFFERENT  VALUES 
OF  B  .  TESTS  COMPUTED  WITH  PENTAGONAL  PYRAMID 


FIGURE  8.15  -  NORMALIZED  BUBBLE  PERIOD  VERSUS 
INVERSE  NORMALIZED  BUBBLE-APEX 
DISTANCE 


FIGURE  8.17  -  MOTION  OF  BUBBLE  CENTER  VERSUS  TIME  FOR  DIFFERENT 
VALUES  OF  6  .  TESTS  CONDUCTED  WITH  2-D  CORNER 
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FIGURE  8.18  -  MOTION  OF  BUBBLE  CENTER  VERSUS  TIME  FOR  DIFFERENT 
VALUES  OF  6  .  TESTS  CONDUCTED  WITH  3-D  CORNER 


FIGURE  8.19  -  MOTION  OF  BUBBLE  CENTER  VERSUS  TIME  FOR  VALUES 
OF  6  .  TSET  CONDUCTED  WITH  PENTAGONAL  PYRAMID 


161*0 


FIGURE  8.20  -  MOTION  OF  BUBBLE  CENTER  VERSUS  TIME  FOR  DIFFERENT 
VALUES  OF  &.  TESTS  CONDUCTED  WITH  CONE 
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FIGURE  8.21  -  GROWTH  AND  COLLAPSE  OF  A  SPARK  GENERATED  BUBBLE 

INSIDE  CONE.  d$  =  4  inches,  AP  =  1.4  x  104Pa,  TIME 

PER  FRAME  =  0.  249  ms. 
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FIGURE  8.22  - 


GROWTH  AND  COLLAPSE  OF  SPARK  GENERATED  BUBBLE 
INSIDE  CONE.  ds  =  H  inches,  ^P  =  1.6  x  104Pa,  TIME 

PER  FRAME  =  0.294  ms. 
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FIGURE  9.1  -  DIAGRAM  OF  EXPERIMENTAL  SET-UP 


FIGURE  9.2  -  PHOTOGRAPH  OF  WATER  TANK  USED  FOR  ACOUSTIC 

DAMPING  TESTS  SHOWING  BOTH  HYDROPHONES  AND 
A  BUBBLE  SCREEN 


PERCENT  SIGNAL  DAMPED 
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FIGURE  9.4  -  CALIBRATED  FREQUNCY  RESPONSE  OF  CERAMIC 
CYLINDER  HYDROPHONE,  (TH-653)  AND  KSP- 
UT-114  HYDROPHONE 
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FIGURE  9.5  -  DAMPING  RATIO  VERSUS  EMISSION  FREQUENCY  AT  CONSTANT 
TRANSMITTING  HYDROPHONE  INPUT  VOLTAGE  COMPARED  WITH 
THE  RECEIVED  SIGNAL  IN  ABSENCE  OF  BUBBLES 


SIGNAL  AMPLITUDE  -  NO  BUBBLES,  VOLTS 


Tracor  Hydronautics 


NO  BUBBLES  -  CONSTANT  OUTPUT 


AIR  INJECTION  PRESSURE:  10  psig 
RECEIVING  HYDROPHONES:  KSP-UT-114 


FREQUENCY,  kHz 


FIGURE  9.6  -  DAMPING  RATIO  VERSUS  EMISSION  FREQUENCY 
VARIABLE  INPUT  VOLTAGE 


FIGURE  9.7  -  PICTURE  OF  BUBBLE  STREAM  GENERATED  WITH 
HYDROPERM  TUBES.  AIR  INJECTION  PRESSURE: 
10  psiq .  WIRE  SIZES:  5,  10,  AND  15  mils. 


PICTURE  OF  BUBBLE  STREAM  GENERATED  WITH 
HYDROPERM  TUBES.  AIR  INJECTION  PRESSURE: 
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FIGURE  9.9  -  BUBBLE  SIZE  DISTRIBUTION  AS  MEASURED  FROM 
PHOTOGRAPHS  AND  DERIVED  FROM  ATTENUATION 
MEASUREMENTS.  HYDROPERM  TUBE  -  10  psig 
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FIGURE  9.11  -  PHOTOGRAPH  OF  BUBBLE  STREAM  GENERATED  WITH 

ALUMINUM  TUBE  WITH  DRILLED  HOLES  0.006  in. 
DIAMETER  -  AIR  INJECTION  PRESSURE:  3  psiq 
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INJECTION  AIR  PRESSURE  =  10  psig 
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FIGURE  9.13  -  BUBBLE  SIZE  DISTRIBUTION  AS  MEASURED  FROM 

PHOTOGRAPHS  AND  DERIVED  FROM  ATTENUATION 
MEASUREMENTS.  ALUMINUM  TUBE  WITH  0.06  In. 
DIAMETER  HOLES 
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FIGURE  9. H  -  DAMPING  RATIO  VERSUS  FREQUENCY  FOR  THREE 
REPEATED  TESTS 


